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ABSTRACT 


A new  hydrodynamical  interpolation  technique  is  developed  and  tested  to  construct  a model  of  global  ocean 
tides  with  the  support  of  empirical  tidal  constants  collected  around  the  world.  The  discrete  tide  model  features  a 
1°  by  1°  spherically  graded  grid  system  in  connection  with  a hydrodynamically  defined  bathymetry  recognizing 
barrier  effects  of  large  boundary  and  bottom  anomalies.  The  Laplace  tidal  equations  are  augmented  by  turbulent 
friction  terms  with  novel  mesh-area  (latitude  and  depth)  dependent  eddy-viscosity  and  bottom-friction  coefficients. 
The  well-known  astronomical  tide-generating  forces  are  modified  by  simplified  versions  of  secondary  effects  due  to 
solid  earth  tides  and  ocean-tidal  loading.  Central  staggered  finite  differences  in  space  and  new  averaged  finite  dif- 
ferences in  time  are  used  to  enhance  decay,  dispersion,  and  stability  characteristics  and  to  facilitate  the  hydro- 
dynamical  interpolation  of  empirical  data.  This  unique  interpolation  is  accomplished  by  a controlled  adjustment  of 
the  bottom-friction  coefficient  and  by  redefining  a more  physical  shoreline  through  a monitored  in-  or  out-flow 
allowed  across  the  artificial  mathematical  coastline. 

Extensive  -computer  experiments  were  conducted  to  study  the  characteristics  of  the  novel  friction  laws  and 
hydrodynamical  interpolation  methods  and  to  determine  by  trial-and-error  optimum  friction,  interpolation,  and 
finite-differencing  parameters.  The  computed  M2-tide  data  along  with  all  (specially  labeled)  empirical  constants  are 
tabulated  in  geographical  map  form  for  four  typical  30°  by  50°  ocean  areas.  A complete  tabulation  and  discussion  of 
the  computed  M2  tide  will  be  published  in  Part  II  of  this  report.  It  is  estimated  that  the  tabulated  tidal  charts  permit 
a prediction  of  the  Mj-tide  elevation  of  the  ocean  surface  over  the  geoidal  rest  level  with  an  accuracy  of  better  than 
5 cm  anywhere  in  the  open  ocean  and  with  somewhat  less  accuracy  near  rough  shorelines.  With  the  forthcoming 
construction  of  the  lesser  S2,  N2,  and  K2;  Kj,  O, , Pj,  and  Q( ; and  Mf,  Mm,  and  Ssa  tidal  constituents,  the  total 
tide-prediction  error  can  be  kept  below  the  10-cm  bound  posed  by  applied  researchers  of  today. 
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1.  INTRODUCTION 

The  restless  sea  with  its  surfing  waves,  swirling  currents,  and  swelling  tides  has  fascinated,  worried,  and  fright- 
ened mankind  from  the  beginning.  The  mighty  erosive  and  devastating  power  of  ocean  tides  flooding  shorelands  and 
seaports  always  troubled  land  developers,  harbor  engineers,  and  seamen.  The  observed  periodic  feature  of  tides  very 
early  led  laymen  tidalists  to  record,  analyze,  and  predict  high  and  low  tides  in  coastal  waters  and  to  warn  and  pre- 
pare people  to  take  appropriate  precautions.  Indeed,  quite  satisfactory  tide  tables  were  constructed  by  simple  rules 
of  thumb.  Nowadays,  in  seashore  areas,  newspapers,  radio  and  television  broadcasts,  and  special  bulletins  publish 
daily,  monthly,  and  longer-range  tide  predictions  of  interest  to  many  users. 

Up  to  recent  years,  practical  interest  in  ocean  tides  was  essentially  confined  to  coastal  waters.  With  the  ad- 
vancement of  science  and  technology,  the  need  for  extremely  accurate  tide  predictions  in  all  the  world  oceans  has 
become  an  urgent  problem.  In  fact,  ocean  tides  represent  fluctuating  loads  on  the  solid  earth,  which  cause  tilting  of 
the  crust  and  disturbances  of  the  earth's  stress  and  gravity  fields.  The  precise  knowledge  of  these  loads  permits 
researchers  to  determine  important  hydrodynamical  parameters  of  the  oceans  and  elastic  parameters  of  the  solid 
earth.  Similarly,  interactions  between  ocean  tides  and  atmosphere,  celestial  bodies,  and  artificial  satellites  and 
missiles  can  be  studied  with  high  accuracy. 

Recently,  the  National  Aeronautics  and  Space  Administration  (NASA)  and  the  Department  of  Defense  (DoD) 
joined  in  the  common  mission  to  map  the  geoid  at  sea  by  satellite  altimetry  up  to  about  10-cm  accuracy.  However, 
the  geoidal  (rest)  sea  level  is  hydrodynamically  disturbed  by  ocean  tides  and  other  currents  due  to  pelagic  density 
variations  and  atmospheric  surface  forces.  Such  hydrodynamical  undulations  of  the  sea  surface  topography  reach 
amplitudes  of  more  than  100  cm  in  open  oceans  and  200  cm  in  coastal  waters.  Thus,  ocean  tides  and  currents  need 
to  be  determined  with  an  accuracy  compatible  with  the  desired  geoid  accuracy  in  order  to  provide  effective  correc- 
tions for  altimeter  measurements. 

In  brief,  contemporary  researchers  of  today  in  geophysics,  geodesy,  oceanography,  meteorology,  astronomy, 
and  space  technology  all  pose  the  same  challenging  question: 

We  need  to  know,  for  any  given  time  and  place  of  the  world  oceans,  the  tidal  elevation  of  the  sea  surface 

over  its  geoidal  level  within  10-cm  accuracy. 

In  the  past  two  decades,  ocean  tidalists  have  devoted  a considerable  effort  to  answer  this  extraordinary  question  by: 

a.  Empirical  methods;  that  is.  by  means  of  vast  numbers  of  tidal  recordings  taken  at  continental,  island,  and 
deep-sea  stations  around  the  world. 

b.  Theoretical  methods;  i.e.,  by  tapping  the  wealth  of  novel  mathematical  tools  designed  by  computer-ori- 
ented numerical  analysts  for  the  integration  of  appropriate  hydrodynamical  equations  governing  tidal  currents  with 
theoretical  and/or  empirical  boundary  and  initial  conditions. 

Considerable  progress  in  the  qualitative  mapping  of  ocean  tides  has  indeed  been  made.  Many  realistic  features 
of  oceanic  tides  have  been  discovered  and  interpreted.  This  rich  and  extensive  work  will  be  briefly  discussed  in 
Section  2,  Reviews  and  Previews.  However,  in  spite  of  the  massive  empirical  and  theoretical  effort,  the  satisfactory 
answer  to  the  question  posed  above  remained  elusive.  The  constructed  tidal  charts  vary  over  large  ocean  areas  from 
investigator  to  investigator,  and  tide  predictions  fall  considerably  short  of  the  desired  accuracy.  Nevertheless,  the 
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published  tidal  maps,  for  instance,  those  by  Pekeris  and  Aecad  (1969).*  Zahel  (1970,  1973,  1975).  and  Estes  (1975, 
1977),  do  suggest  that  the  desired  mapping  of  ocean  tides  might  succeed  with  somewhat  more  refined  models.  A 
similar  situation  was  recognized  (Doiis  et  al„  1974;  and  Reid  ei  al„  1975)  by  numerical  oceanographers  and  mete- 
orologists modeling  even  more  involved  general  ocean  currents  and  atmospheric  circulations.  Obviously,  the  success- 
ful study  of  such  turbulent  fluid  motions  depends  most  of  all  on  the  realistic  representation  of  the  dissipative  forces 
acting  in  the  fluid  interior  and  along  its  boundary.  The  special  field  of  ocean  tides,  which  are  distinguished  by  their 
simplifying  periodic  property  and  by  their  widely  observed  boundary  values,  may  be  considered  a proving  ground  for 
other  models  of  turbulent  fluid  motions. 

In  the  following  treatise,  a new  combined  hydrodynamical  and  empirical  technique  will  be  introduced  and 
tested  to  construct,  at  first,  the  principal  semidiurnal  moon  (M2)  tide.  This  procedure  represents  a modified  version 
of  the  well-known  hydrodynamical  numerical  method  developed  and  tested  by  Hansen  (1966)  and  his  coworkers 
Friedrich  (1966),  Brettschneider  (1967),  Trepka  (1967),  and  Zahel  (1970,  1973.  1975),  as  well  as  by  Estes  ( 1 975 , 
1977),  Among  other  modifications,  the  most  significant  improvements  of  the  results  were  attributed  to: 

a.  A unique  “hydrodynamical  interpolation"  of  more  than  2 000  empirical  tidal  data  into  the  original  purely 
hydrodynamical  tide  model.  This  interpolation  is  accomplished  by  a controlled  local  adjustment  of  the  bottom- 
friction  coefficient  and  by  allowing  a monitored  in-  or  out-flow  across  the  mathematical  ocean  boundary  and,  thus, 
redefining  implicitly  a more  physical  boundary. 

b.  A graded  1°  by  1°  grid  system,  in  connection  with  a hydrodynamically  defined  bathymetry  recognizing 
the  important  barrier  effects  of  narrow  ocean  ridges. 

c.  New  averaged  finite  differences  in  time  enhancing  decay,  dispersion,  and  stability  characteristics  and 
facilitating  the  implicit  hydrodynamical  interpolation  of  empirical  data  (a). 

d.  A novel  mesh-area  (mesh-size  and  depth)  dependent  eddy  viscosity  specifying  a more  realistic  eddy  dis- 
sipation. 

The  complete  results  for  the  constructed  M2  tide  will  be  published  in  table  and  map  form  in  Part  II  of  this 
report  (Schwiderski,  1978  b).  Similar  charts  for  the  S2,  N2,  and  K2;  K, , O, . P, , and  0, ; and.  possibly.  Mf,  Mm, 
and  Ssa  tides  (Table  I)  will  be  constructed  and  published  as  additional  parts  of  this  report.  A separate  tabulation 
(Schwiderski.  1978  a)  of  the  new  hydrodynamically  defined  ocean  bathmetry  is  in  preparation.  All  tidal  and  bathy- 
metry data  will  be  available  in  tape  form  at  NSWC,  Dahlgren,  Virginia. 


*A  list  of  references  is  included  at  the  end  of  this  report. 
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Tabic  I . Constants  of  Major  Tidal  Modes 


Tidal  Mode 

*a(m) 

ob(10'4/sec) 

Xf(deg) 

Semidiurnal  Species 

M,  = Principal  Lunar 

0.242  334 

1.405  19 

2h0  - 2s0 

S2  = Principal  Solar 

0.113  033 

1 .454  44 

0 

N2  = Elliptical  Lunar 

0.046  398 

1.378  80 

2h0  - 3sq  + Po 

K-,  = Declination  Luni-Solar 

0.030  704 

1 .458  42 

2h0 

Diurnal  Species 

= Declination  Luni-Solar 

0.141  565 

0.729  21 

h0  ♦ 90 

O,  = Principal  Lunar 

0.100  514 

0.675  98 

h0  - 2s0  - 90 

Pj  = Principal  Solar 

0.046  843 

0.725  23 

-h0  - 90 

Qj  = Elliptical  Lunar 

0.019  256 

0.649  59 

ho  - 3so  + Po  - 

Long-Period  Species 

Mf  = Fortnightly  Lunar 

0.041  742 

0.053  234 

2*0 

Mm  = Monthly  Lunar 

0.022  026 

0.026  392 

0 
Q. 

1 

O 

t/i 

Ssa  = Semiannual  Solar 

0.019  446 

0.003  982 

2h0 

aA.-  = amplitude  of  the  partial  tide. 
bo  = frequency  of  the  partial  tide. 

CX  = astronomical  argument  of  the  partial  tide. 

^(hg,  Sq.Pq)  - mean  longitudes  of  sun.  moon,  and  lunar  perigee  at  Greenwich  midnight; 
hQ  = 279.696  68  + 36  000.  768  930  485  T + 3.03-10'4T2. 
sQ  = 270.434  358  + 481  267.883  141  37T  - 0.001  133TJ  + 1.9*10-‘T3. 

P0  = 334.329  653  + 4 069.034  032  957  5T  - 0.010  325T2  - 1.2-10‘5T3, 
where 

T = [27  392.500  528  + 1 .000  000  035  6D) /36  525, 

D = d + 365(y  - 1975)  + Int  |(y  - 1975)/4[. 
d = day  number  of  year  (d  = 1 for  January  1 ). 
y > 1975  = year  number, 
and 

Int  |x|  = integral  part  of  x. 
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2.  REVIEW  AND  PREVIEW 


Very  early,  scientists  linked  the  periodic  phenomenon  of  ocean  tides  to  the  apparent  motions  of  the  moon 
and  sun.  The  first  physical  explana  i of  tides  was  given  by  Newton  (1687),  who  applied  his  newly  discovered  theo- 
ry of  gravitation  and  introduced  the  still  important  “equilibrium  theory”  of  tides  (Sect.  3.D):  under  the  attraction 
of  the  moon  and  the  sun,  the  sea  surface  is  assumed  to  attain  instantaneously  the  shape  of  a level  surface  in  hydro- 
static equilibrium.  This  theory  was  perfected  after  Newton  in  1738  (Cantor,  1901)  in  a contest  conducted  by  the 
Paris  Academy  of  Sciences.  Among  the  leading  contestants  for  the  best  mathematical  and  physical  description  of 
ocean  tides  were  such  distinguished  scientists  as  D.  Bernoulli,  Euler,  and  Maclaurin.  The  hydrostatic  equilibrium 
theory  explained  in  fact  the  periodicity  and  some  other  observed  simple  features  of  ocean  tides,  but  failed  entirely 
to  describe  all  hydrodynamical  properties  that  were  clearly  observed. 

Although  Newton  recognized  the  distinct  hydrodynamical  nature  of  ocean  tides,  it  was  Laplace  (1775)  who 
formulated  the  first  hydrodynamical  equations  of  oceanic  tidal  motions  (Sect.  3.F).  The  Laplace  tidal  equations 
(LTEs)  consider  an  inviscid  and  incompressible  fluid  subject  to  kinetic,  potential,  and  Coriolis  forces  generated  by 
the  primary  astronomical  tide-producing  potential,  which  is  directly  proportional  to  the  equilibrium  tide.  Numerous 
solutions  of  these  equations  were  constructed  by  analytic  and  semi-analytic  methods  for  idealized  ocean  basins, 
estuaries,  and  channels.  Realistic  results  were  sought  by  many  well-known  researchers,  such  as  Laplace  (1775), 
Airy  (1842),  Ferrel  (1874),  Thomson  (Lord  Kelvin.  1879),  Hough  (1897),  Poincare  (1910),  Lamb  (19L),  Proud- 
man  (1915),  Proudman  and  Doodson  (1927),  Platzman  (1971 ; 1972  a,  b;  1974),  and  Miles  ( 1 974). 

Most  of  this  rich  and  elucidating  work  has  been  collected  and  analyzed  in  books  and  reviews  (see  References). 
The  conclusions  clarified  many  realistic  properties  of  ocean  tides  including  the  existence  of  so-called  amphidromic 
points  around  which  oceanic  tidal  waves  rotate  with  amazing  wave  speeds  (Sec.  6.B).  Nevertheless,  the  accurate  and 
detailed  prediction  of  ocean  tides  remained  unsolved.  More  realistic  physical  and  geometrical  models  and  advanced 
mathematical  techniques  were  obviously  needed  to  compute  the  desired  ocean  tides. 

The  theoretical  mouels  as  well  as  the  empirical  predictions  of  ocean  tides  were  greatly  aided  by  Thomson 
(Lord  Kelvin,  1868),  who,  following  an  earlier  suggestion  by  Laplace,  introduced  the  method  of  harmonic  analysis 
into  tidal  studies.  The  astronomical  tide-generating  potential  (equilibrium  tide)  is  expanded  into  an  almost  periodic 
(with  a nonharmonic  frequency  spectrum)  series  of  harmonic  tidal  components  (Sect.  3.D).  Later,  G.  H.  Darwin 
(1883)  improved  the  decomposition  up  to  39  terms,  and  Doodson  (1921)  carried  it  on  to  400  frequency  modes  by 
introducing  his  ingenious  Doodson  numbers.  A novel  nonharmonic  and,  subsequently,  harmonic  decomposition  of 
the  equilibrium  tide  was  derived  by  Cartwright  and  Tayler  ( 1970). 

The  harmonic  expansion  of  the  tide-generating  potential  is  evidently  of  fundamental  significance  for  the 
mathematical  analysis  and  prediction  ol  ocean  tides.  If  one  assumes  that  ocean  tides  arc  governed  by  linear  or  almost 
linear  equations  of  motion  (e.g.,  LTEs,  Sect.  3.F),  then  the  oceanic  tide  is  also  representable  as  a linear  superposition 
of  harmonic  modes.  Every  harmonic  component  of  the  equilibrium  tide  generates,  through  the  ocean’s  response,  a 
similar  component  of  the  oceanic  tide  of  identical  frequency  differing  only  in  two  time-independent  constants:  the 
amplitudes  and  phases.  Once  these  harmonic  constants  are  determined  for  a given  geographical  location  either  by 
observation  or  by  theoretical  computations,  then  the  time-dependent  tide  can  be  predicted  by  simple  linear  super- 
position. Moreover,  every  oceanic  tidal  mode  is  hydrodynamically  independent  of  all  others  and  can  be  constructed 
without  any  knowledge  of  the  others. 

However,  some  nonlinear  oceanic  responses  in  coastal  areas,  where  the  sea  surface  varies  rapidly,  may  be 
necessary  in  a realistic  tide  model  (Sect.  5.F).  There,  one  must  remember  that  any  nonlinearity  generates  higher 
harmonic  frequencies  and  causes  involved  interactions  between  different  tidal  modes  and  also  other  ocean  currents. 
In  this  connection,  one  may  compare  the  introductory  remarks  to  the  British  Admiralty  (1977)  tide  tables. 
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Over  many  years,  massive  tidal  records  (some  longer  than  half  a century  with  more  than  10  million  hourly 
readings)  have  been  made  with  about  10  000  tide  gauges  placed  around  the  world  at  continental  and  island  stations. 
Most  of  these  data  have  been  harmonically  analyzed  for  many  ocean  tide  components.  The  extracted  harmonic  con- 
stants have  been  tabulated  and  are  now  available  at,  for  instance,  the  National  Ocean  Survey  (1942),  the  Interna- 
tional Hydrographic  Bureau  (1966),  and  the  British  Admiralty  (1977). 

In  recent  years.  Eyrie  (1968),  Filloux  (1968),  and  Snodgrass  (1968)  developed  sensitive  deep-sea  tide  (pres- 
sure) gauges  that  can  be  placed  almost  anywhere  on  the  ocean  floor  to  record  tidal  variations  for  up  to  one  year. 
Munk  and  Cartwright  (1966)  introduced  the  so-called  “ocean  response  method”  in  order  to  analyze  relatively  short- 
time  records  with  sufficient  accuracy.  Though  this  method  is  basically  nonharmonic,  it  too  yields  as  its  end  product 
the  most  desired  harmonic  constants.  The  local  time  series  of  the  ocean  tide  is  considered  as  a convolution  of  the 
tide-generating  potential  and  a transfer  function  (complex  response  weights)  which  is,  a priori,  not  known  and, 
hence,  must  be  determined  in  an  optimum  sense  (Cartwright,  1968,  1969;  and  Cartwright  et  al.,  1969).  By  using 
a..alyzed  coastal  reference  stations,  satisfactory  results  are  being  claimed  for  rather  short  recording  times  of  two 
weeks  and  even  one  week.  Deep-sea  measurements  have  been  published,  for  instance,  by  Munk  et  al.  (1970),  Irish 
et  al.  (1971),  Nowroozi  (1972),  Mofjeld  (1975),  Pearson  (1975  a,  b),  and  Zettler  et  al.  (1975). 

By  simple  inspection  of  numerous  observed  ocean-tide  data  and  experienced  intuition,  cotidal  (equi-phase) 
and  some  corange  (equi-amplitude)  maps  have  been  constructed  for  marginal  seas  and  partial  or  worldwide  oceans 
by  Whewell  (1833,  1848),  Berghaus  (1845),  Harris  (1904),  Sterneck  (1920, 1921),  Prufer  (1939»,  Dietach  (1944  a, 
b),  Villain  (1952),  Bogdanov  (1961  a,  b),  Luther  and  Wunsch  (1974),  and  others.  All  of  these  map;  give  a qualitative 
overview  of  the  phenomena  of  ocean  tides. 

A more  systematic  empirical  method  to  map  large-scale  ocean  tides  was  proposed  and  demonstrated  by  Kuo 
et  al.  (1970  a,  b),  .Jachens  and  Kuo  (1973),  and  Kuo  and  Jachens  (1977).  In  their  inversion  method,  modern  optimi- 
zation techniques  are  utilized  to  fit  low-degree  polynomials  simultaneously  to  observed  ocean-tide  data  and  to 
measured  oceanic  tidal  disturbances  of  the  solid-earth  gravity,  which  are  correlated  through  well-known  interaction 
equations.  While  this  approach  is  definitely  more  reliable  than  purely  intuitive  attempts,  it  appears  to  be  constrained 
by  rather  low-degree  polynomials  which  can  hardly  accommodate  the  various  tidal  undulations  known  from  obser- 
vations (Sect.  6.B).  The  method  also  offers  no  physical  insight  into  the  tidal  phenomena  and  yields  no  information 
on  tidal  currents,  which  is  needed  in  certain  applications.  Aside  from  those  shortcomings,  the  method  seems  to  yield 
promising  results  in  limited  ocean  regions. 

With  the  advent  of  electronic  computers,  hydrodynamical  numerical  methods  dominated  the  theoretical  scene 
to  map  worldwide  ocean  tides.  The  modern,  large-scale  numerical  investigation  of  ocean  tides  was  pioneered  by 
Hansen  (1948,  1949,  1966),  who  realized  that  the  complexity  of  natural  ocean  basins  renders  any  analytic  treatment 
intractable.  Hansen  was  soon  followed  by  a long  line  of  researchers,  such  as  Gohin  (1961),  Accad  and  Pekeris 
(1963),  Ueno  (1963,  1964),  Bogdanov  et  al.  (1964),  Bogdanov  and  Magarek  (1967,  1969),  Brettschneider  (1967), 
Tiron  (1967),  Trepka  (1967),  Pekeris  and  Accad  (1969),  Marchuk  et  al.  (1969),  Zahel(1970,  1973,  1975),  Hender- 
shott  (1972,  1975),  Gordeyev  et  al.  (1973),  Marchuk  et  al.  (1973),  and  Estes(1975,  1977).  Their  probing  computer 
experiments  of  steadily  increasing  sophistication  were  paralleled  and  supported  by  an  even  longer  list  of  numerical 
analysts  modeling  other  fluid  motions,  in  particular  closely  related  general  ocean  currents  and  atmospheric  circu- 
lations. Among  the  many  investigators,  only  very  few  may  be  named  in  this  connection  because  of  their  important 
contributions  with  direct  relevance  to  ocean  tides:  Bryan  (1963),  Sma.’orinsky  (1963),  Friedrich  (1966,  1970), 
Crowly  (1968,  1970),  Leith  (1968),  Cox  (1970),  O’Brien  (1971),  and  Holland  and  Hirschman  (1972),  as  well  as 
Dooset  al.  (1974),  and  Reid  et  al.  (1975). 


Hansen  (1948,  1949)  began  his  large-scale  hydrodynamical  numerical  investigations  of  ocean  tides  with  the 
classical  LTEs  (Sect.  3.F).  These  equations  are  derived  from  the  Euler-Lagrange  equations  of  inviscid  and  incompres- 
sible fluid  motions  by  integration  over  the  ocean  depth,  invoking  the  hydrostatic-pressure  assumption  and  neglecting 
all  nonlinear  effects  and  any  motion  in  depth  (Sect.  3.E).  Most  of  those  simplifying  conditions  could  be  justified  by 
long-wave  theory.  The  basic  LTEs  are  the  starting  point  for  all  ocean  tide  investigators.  The  quality  of  their  work 
depends  essentially  on  their  more  or  less  realistic  modeling  of  the  following  hydrodynamical  and  mathematical 
components  characterizing  true  tidal  motions. 

(a)  Lateral  Eddy  Dissipation.  To  improve  his  results,  Hansen  (1966)  recognized  the  fundamentally  turbulent 
nature  of  tidal  currents  and  augmented  the  LTEs  by  lateral  eddy -dissipation  terms  (Sect.  3.C).  Such  interior  friction 
forces  were  first  introduced  into  fluid  mechanics  by  Boussinesq  (1877),  who  replaced  the  unknown  turbulent 
Reynolds  stress  tensor  by  the  known  laminar  stress  tensor  with  a constant  eddy  viscosity  at  one’s  disposal.  In 
oceanography,  this  simple  concept  was  brought  into  practice  by  Ekman  (1902)  and  was  later  effectively  used,  for 
instance,  by  Munk  (1950)  and  Bryan  (1963). 

Ignoring  the  successes  of  Boussinesq’s  unique  substitution,  which  serves  to  control  the  supercritical  instability 
of  turbulent  flow  in  place  of  the  Reynolds  stresses,  many  researchers  still  regard  this  approximation  as  controversial, 
artificial,  fictitious,  and  phony.  Apparently,  this  difficulty  originates  from  the  observation  that  the  “constant”  eddy 
viscosity  varies  enormously  in  form  and  magnitude  from  investigator  to  investigator.  While  the  kinematic  molecular 
viscosity  of  water  is  of  the  order  10"2  cm2/sec,  oceanic  eddy  viscosities  of  the  order  10°  to  1011  cm2/sec  have  been 
quoted  (Cox,  1970;  Zahel,  1970).  Nevertneless,  the  huge  variation  of  the  eddy  viscosity  is  physically  plausible,  and 
some  empirical  evidence  for  the  lower  value  has  been  discussed  by  Munk  (1966). 

As  was  first  realized  by  Prandtl  (1925),  'ne  eddy  viscosity  is  a characteristic  quantity  of  the  turbulent  motion 
under  investigation  and  not  just  a fluid  cc-.stant  such  as  the  molecular  viscosity.  In  fact,  a prescribed  eddy  viscosity 
actually  specifies  the  “macroscopical  (averaged)  features  of  the  “microscopically  undefined”  turbulent  flow. 
Consequently,  the  quality  of  the  modeling  of  any  turbulent  motion  is  directly  linked  to  the  investigator's  realistic 
choice  of  the  eddy  viscosity. 

The  physically  sound  meaning  of  Boussinesq’s  approximation  of  the  Reynolds  stresses  was  first  illustrated  by 
Prandtl’s  (1925)  celebrated  "momentum  austausch  (exchange,  mixing)  theory”  of  turbulent  flow.  With  the  help  of 
constant  “mixing-length  parameters"  (see.  e.g.,  Schlichting,  1968),  the  eddy  viscosity  (momentum  austausch  coef- 
ficient) is  found  to  be  dependent  on  the  averaged  flow  velocity;  i.e.,  the  eddy  dissipation  appears  as  a quadratic 
phenomenon.  In  meteorology  and  oceanography,  velocity-dependent  eddy  viscosities  were  derived  on  physical 
grounds  by  Smagorinsky  (1963),  Leith  (1968),  Crowley  (1968,  1970),  O’Brien  (1971),  and  others. 

Considering  the  complexity  of  the  world  oceans,  it  is  attractive  to  retain  the  simplifying  linearity  of  the  tidal 
equations.  Since  averaged  tidal  velocities  can  be  expected  to  be  generally  slow,  it  appears  more  consistent  with  the 
overall  Stokes  slow-motion  assumption  underlying  the  basic  LTEs  to  employ  a velocity-independent  eddy  viscosity 
(compare  bottom  friction  below).  Also,  the  refined  studies  of  Zahel  (1970.  1973,  1975),  Marchuk  et  al.  (1973),  and 
Estes  (1975,  1977)  do  not  exhibit  any  need  for  nonlinear  eddy -dissipation  terms.  However,  the  author's  frustrating 
early  experiments  with  a worldwide,  constant  eddy  viscosity  in  an  ocean  of  depths  varying  from  10  to  7 000  m 
(Sect.  5.B)  clearly  indicated  that  the  lateral  eddy  viscosity  must  depend  linearly  on  the  ocean  depth  (Sect.  3.C). 
Eddy  viscosities  depending  on  1/2  and  3/2  powers  of  the  depth  were  also  tried,  but  they  failed  to  yield  equally  good 
results. 


A renewed  look  at  the  phenomenon  of  turbulent  motion  led  to  the  introduction  of  a lateral  eddy  viscosity 
that  depends  linearly  on  a horizontal  (cell  or  mesh  size)  and  a vertical  (cell  size  or  depth)  mixing  length  (Eq.  6). 


This  physically  plausible  law  explains  the  enormous  variations  of  eddy  viscosities  mentioned  above.  The  fact  that 
the  eddy  viscosity  must  depend  in  some  way  on  the  mesh  size  (resolution)  in  discrete  models  of  ocean  currents  was 
earlier  noticed  by  Cox  (1970).  Friedrich  (1970).  Holland  and  Hirschntan  (1972).  and  Zahel  (1975).  It  may  be 
mentioned,  that  the  author's  extensive  computer  experiments  with  the  novel  law  of  eddy  dissipation  resulted  in 
rather  drastic  improvements  of  the  ocean-tide  model. 

(bl  Bottom-Friction  Law.  In  addition  to  eddy  dissipation.  Hansen  (1966)  supplemented  the  LTEs  also  by 
quadratic  bottom-friction  terms.  In  fluid  mechanics,  the  physical  significance  of  the  quadratic  law  of  wall  friction 
was  first  pointed  out  by  Boussinesq  (1896).  Later.  G.  1.  Taylor  (1918)  applied  it  advantageously  in  his  study  of 
tidal  currents  in  the  Irish  Sea. 

In  order  to  retain  the  linearity  of  the  tidal  equations  for  their  adopted  time-independent  numerical  procedure, 
Pekeris  and  Accad  (1969)  used  the  linear  law  of  bottom  friction  with  a coefficient  depending  inversely  on  the  ocean 
depth.  These  authors  also  called  attention  to  the  earlier  work  by  Grace  (1931).  who  applied  both  the  quadratic  and 
the  linear  laws  of  bottom  friction  to  the  problem  of  tides  in  the  Gulf  of  Suez.  Grace  experienced  a slight-preference 
for  the  linear  law.  Similarly  to  eddy  dissipation  (see  above),  it  appears  that  the  linear  law  of  bottom  friction  is  more 
consistent  with  all  other  assumed  linearizations  of  the  equations  of  slow  (Stokes-like)  averaged  tidal  motion  (see, 
e.g.,  Schlichting  (1968)).  Furthermore,  the  convergence  characterisitcs  of  the  time-stepping  computations  carried 
out  and  displayed  by  Estes  (1975)  with  the  quadratic  law  of  bottom  friction  fail  to  exhibit  any  nonlinear  symptoms. 

Other  bottom-friction  laws  were  considered  by  Johns  (1966)  and  McGregor  (1972),  as  well  as  by  Kagan 
(1971  a.  b).  In  Section  3.B,  the  linear  law  of  bottom  friction  will  be  selected  as  preferable  to  the  quadratic  law  at 
least  over  all  open-ocean  areas.  In  contrast  to  Pekeris  and  Accad  (1969),  the  bottom-friction  coefficient  was  found 
more  realistic  without  any  dependence  on  the  ocean  depth,  which  varies  in  the  present  model  from  10  to  7 000  m. 
Analogous  to  the  eddy  viscosity  (see  (a)),  the  bottom-friction  coefficient  (as  a vertical  eddy  viscosity)  is  assumed  on 
physical  arguments  to  depend  on  the  bottom  cell  area;  that  is.  on  the  two  horizontal  mesh  sizes  of  the  discrete  tidal 
model  considered  (Eq.  4).  Certain  nonlinear  bottom-friction  effects  will  be  introduced  indirectly  along  coastal 
boundaries  by  the  novel  technique  of  hydrodynamical  interpolation  of  empirical  tide  data  described  below  and  in 
Section  5.F. 

(c ) Primary  and  Secondary  Tide-Generating  Forces.  Most  numerical  tidalists  use  only  the  moon's  or  sun  s 
tide-generating  potential  as  the  primary  driving  force  of  ocean  tides.  For  obvious  mathematical  advantages  (see 
above),  only  one  harmonic  (mostly  the  largest  M,)  component  is  investigated  (Sect.  3.D).  Recently  recognized 
significant  interactions  between  the  solid  earth  and  ocean  tides  motivated  Farrell  ( 1972  a.  b;  1973)  to  suggest  that 
accurate  models  of  ocean  tides  should  include  tidal  deformations  of  the  solid  earth,  tidal  loading  of  the  ocean,  and 
their  associated  perturbations  of  the  primary  tide-generating  potential.  Farrell  went  on  to  derive  the  mathematical 
expressions  of  those  secondary  tide-modifying  forces,  which  should  be  incorporated  in  the  forcing  terms  of  the 
tidal  equations. 

Because  of  the  elastic  properties  of  the  solid  earth,  it  is  sufficient  for  ocean-tide  models  to  include  only 
second-order  Love  number  approximations  (Sect.  3.D)  of  the  terrestrial  tide  and  its  associated  gravity  disturbance. 
This  simple  static  approximation  sets  both  solid-earth  tidal  effects  directly  proportional  to  the  equilibrium  tide. 
Hence,  any  linear  (or  almost  linear)  ocean-tide  model  with  homogeneous  boundary  data  but  without  those  second- 
dary  driving  forces  can  be  corrected,  a posteriori,  by  applying  a uniform  factor  to  all  amplitudes.  No  correction 
of  the  phases  is  required. 

In  order  to  express  the  responses  of  the  solid  earth  and  of  the  gravitational  field  to  the  ocean's  tidal  load, 
Farrell  derived  a Green's  function  expanded  in  spherical  harmonics.  This  representation  changes  the  character  of 
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the  equations  of  oceanic  tidal  motion  from  its  differential  form  to  a much  more  involved  integro-differential  form, 
which  must  be  solved  by  some  iteration  process.  Farrell’s  suggestion  was  picked  up  by  Hendershott  (1972,  1975), 
who  attributed  rather  drastic  effects  to  oceanic  tidal  loading  although  his  numerical  iteration  procedure  failed  to 
converge.  The  same  secondary  forcing  functions  were  also  used  by  Estes  ( 1 977)  in  connection  with  Hansen's  hydro- 
dynamical  numerical  technique  (see  (g)  and  Sect.  5.D)  without  detecting  any  divergence  difficulties.  Estes'  results 
did  not  confirm  the  drastic  effects  of  oceanic  tidal  loading  as  claimed  by  Hendershott.  The  computed  effects  on 
amplitudes  and  phases  are  in  magnitude  considerably  below  those  modifications  that  are  needed  and  that  could 
be  achieved  by  varying  eddy-dissipation  and  bottom-friction  parameters  (Sect's.  3.B  and  3.C).  Nevertheless,  Estes’ 
important  computations  definitely  supported  Farrell’s  original  contention  that  an  accurate  description  of  ocean 
tides  should  include  ocean  loading  effects. 

In  a critical  discussion  of  Farrell’s  representation  of  oceanic  tidal  loading  effects,  Pekeris  (in  a presentation 
given  in  1977  at  the  Office  of  Naval  Research)  noticed  that  the  series  expansion  of  the  Green's  function  is  only 
conditionally  (not  absolutely)  convergent  and  possesses  only  a generalized  derivative  of  Dirac's  6-function  type. 
To  avoid  all  complications  of  this  representation,  Pekeris  suggested  a simple  static  treatment  of  all  effects  due  to 
oceanic  tidal  loading.  Analogous  to  Newton’s  equilibrium  theory  and  to  the  static  approximation  of  the  solid-earth 
tidal  effects  (see  above),  the  solid  earth  and  the  gravity  field  are  assumed  to  respond  instantaneously  to  the  ocean's 
tidal  load;  i.e.,  the  responses  are  directly  proportional  to  the  ocean  tide. 

In  Section  3.D.  Pekeris’  straightforward  solution  of  the  rather  involved  problem  will  be  incorporated  in  the 
present  model.  The  consistent,  approximate  treatment  of  both  effects  due  to  the  solid-earth  tide  and  to  the  ocean’s 
tidal  load  can  be  considered  as  quite  satisfactory,  particularly  for  the  present  ocean  tide  model,  which  will  be  sup- 
ported by  numerous  empirical  tide  data  that  naturally  include  all  solid-earth  and  ocean  responses.  In  agreement  with 
Estes'  results,  the  present  computations  displayed  no  startling  changes  traceable  to  oceanic  tidal  loading. 

(d ) Boundary  and  Initial  Values,  In  order  to  specify  a unique  integral  of  the  LTEs,  Hansen  (1948,  1949) 
supported  his  tide  model  with  empirical  tide  data  as  boundary  values  known  around  the  ocean  basin  considered. 
This  boundary  condition  is  easily  implemented  numerically  by  combining  all  equations  of  oceanic  tidal  motion  into 
a single  elliptic,  second-order  differential  equation  for  the  complex  tidal  amplitude,  provided  a linear  law  of  bottom 
friction  and  no  eddy-dissipation  terms  are  considered.  Although  this  attractive  elliptic  boundary-value  problem  has 
been  studied  by  several  investigators,  it  must  be  emphasized  that  it  leaves  the  fluid  flow  across  the  ocean  boundaries 
without  any  control.  Such  an  unchecked  violation  of  the  stiff  condition  of  conservation  of  mass  is  obviously  some- 
what hazardous.  For  example,  it  is  clearly  one  major  cause  of  the  divergence  problems  encountered  by  Hendershott 
(1975;compare  also  Farrell,  1972  b). 

Instead  of  empirical  tidal  data.  Accad  and  Pekeris  (1963)  and  other  authors  used  the  purely  mathematical 
boundary  condition  of  no-flow  across  the  ocean  shorelines.  When  higher-order  eddy-dissipation  terms  are  included, 
two  sets  of  boundary  data  are  required  to  specify  a solution.  Hansen  (1966),  Zahel  (1970,  1972,  1975),  Estes 
(1975.  1977),  and  others  used  the  mathematical  boundary  conditions  of  no-fiow  across  and  free-slip  along  the 
ocean  boundaries.  Some  authors  (e.g.,  Marchuk  et  al.,  1969)  replaced  the  free-slip  condition  by  the  no-slip ’condi- 
tion, which  is  appropriately  used  to  specify  laminar  flows.  However,  in  turbulent  flows  with  very  thin  boundary 
layers  and  large  eddy  viscosities,  the  free-slip  condition  is  preferred.  Moreover,  this  condition  is  consistent  with  the 
single-layer-ocean  assumption  of  almost  depth-independent  horizontal  velocities.  The  author’s  computer  experiments 
with  both  the  free-slip  and  the  no-slip  conditions  exhibited  a slight  preference  for  the  first  one. 

Some  authors  used  a mixture  of  either  empirical  or  mathematical  boundary  values.  For  example,  Tiron  (1967) 
used  observed  tidal  elevations  wherever  available  and  the  condition  of  o-flow  across  the  coastline  everywhere  else. 
In  any  case,  whenever  empirical  tide  boundary  values  were  incorporated,  all  eddy-dissipation  terms  had  to  be  ne- 
glected and  the  condition  of  conservation  of  mass  had  to  be  violated  without  any  possible  control. 
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In  the  present  tide  model  (Sect.  5.E),  the  mathematical  condition  of  no-flow  across  and  free-slip  along  the 
ocean  boundaries  will  be  strictly  maintained  wherever  no  empirical  tide  data  are  known  (mainly  at  the  coast  of 
Antarctica  and  the  Arctic  Sea).  To  improve  the  quality  of  the  tide  model,  known  empirical  tidal  elevations  will  be 
incorporated  at  some  2 000  shore  points  by  a unique  “hydrodynamica!  interpolation”  technique  without  neglecting 
important  eddy-dissipation  terms  and  permitting  only  a "monitored  small"  break  of  the  no-flow  boundary  condi- 
tion. 


First  attempts  to  use  the  empirical  tide  data  in  place  of  any  mass-conserving  condition  failed  to  produce  realis- 
tic velocity  distributions;  i.e..  unacceptably  large  periodic  mass  flows  in  and  out  of  the  oceans  were  encountered. 
A second  attempt  to  keep  the  empirical  tide  data,  while  enforcing  strictly  the  no-flow  boundary  condition,  by 
adjusting  the  bottom-friction  coefficient  (see  (b)  and  Sect.  5.F)  point-wise  near  the  shore  to  modify  the  velocity 
components,  also  failed  to  yield  realistic  results. 

After  some  evaluation  of  this  disappointing  situation,  it  was  felt  that  quite  satisfactory  results  could  be  ob- 
tained by  a certain  measured  compromise.  The  pointwise  adjustment  of  the  bottom-friction  coefficient  and  the 
resulting  modification  of  the  velocity  field  had  to  be  controlled  by  two  uniform  parameters  (Sect.  5.F),  which 
remained  to  be  determined  by  trial-and-error  computations  to  yield  the  best  results  while  satisfying  the  mass-conser- 
vation condition  as  closely  as  possible.  The  remaining  smaller  violations  of  the  no-flow  boundary  condition  could 
then  be  explained,  for  instance,  by  the  relatively  large  geometric  inaccuracies  of  the  shoreline  of  the  model  (e.g„ 
shape,  si/e.  and  depth  of  boundary  mesh  cells  (see  (0.  Sect.  5.F,  and  Figure  1)).  Naturally,  only  monitored  strictly 
periodic  in-  and  out-flows  could  be  permitted  over  any  boundary  segment,  so  that  the  total  water  mass  of  the 
ocean  remains  conserved.  The  novel  technique  thus  developed  was  found  to  bring  about  the  most  spectacular  im- 
provements of  the  tide  model  desired. 

When  the  complete  time-dependent  problem  is  considered,  its  integration  must  be  started  from  some  meaning- 
ful. initial  oceanic  state,  say.  at  rest.  The  computational  process  can  be  interrupted  at  any  time  for  qualitative  inspec- 
tion and  then  restarted  with  or  without  parameter  changes.  The  integration  can  be  stopped  when  the  time-periodic 
state  imposed  by  the  harmonic  forcing  function  is  reached  with  satisfactory  accuracy.  It  may  be  pointed  out  that 
the  time-dependent  integration  procedure  was  chosen  in  Section  5.D  as  preferable  because  its  stability  analysis 
reveals  important  dispersion,  decay,  and  stability  characteristics  that  determine  the  quality  of  the  tide  model. 
Moreover,  the  experience  gained  with  tidal  time-periodic  currents  can  be  directly  utilized  for  the  modeling  of  more 
involved  general  ocean  currents  and  atmospheric  circulations. 

(e)  Polar  Singularities  and  lie  Caps.  In  spherical  coordinates,  the  equations  of  oceanic  tidal  motion  (COTFs. 
Sect.  3.F I become  singular  at  the  North  Pole.  Most  authors  avoid  this  difficulty  by  terminating  the  ocean  basin  at 
some  fictitious  northern  boundary.  In  order  to  check  his  numerical  method  against  an  analytic  solution  in  a circular 
arctic  ocean  basin.  Zahcl  ( 1970)  derived  an  analytic  integral  of  the  equations  of  motion  without  eddy-dissipation 
terms.  In  Section  4.  a unique  second-order  analytic  solution  of  the  complete  equations  of  tidal  motion  will  be  con- 
structed. This  special  integral  will  be  matched  with  the  numerical  solution  which  starts  at  4°  eolatitude.  It  may  be 
pointed  out  that  the  use  of  the  analytic  North  Pole  solution  allows  for  a more  uniform  grading  (sec  (f)  and  Sect. 
5. A)  of  the  1°  by  1°  grid  system  needed  for  the  definition  of  the  finite  difference  analog.  As  will  be  shown,  a less 
uniform  grid  system  requires  considerably  more  computer  time. 

In  all  tide  models,  it  is  tacitly  assumed  that  the  polar  ice  covers  have  very  little  effect  on  oceanic  tidal  currents. 
This  assumption  is  probably  justified  within  the  desired  accuracy.  It  is  certainly  true  at  least  near  the  North  Pole, 
where,  according  to  the  analytic  solution  (Sect.  4).  all  semidiurnal  and  diurnal  tides  vanish  with  the  second  or  first 
power  ol  the  distance  from  the  Pole.  Furthermore,  the  present  tide  model  is  aided  by  numerous  empirical  tide  data 
collected  at  coastal  stations  of  the  Arctic  Ocean  and  by  some  data  obtained  at  Antarctica,  which  naturally  contain 
ice-cap  effects. 
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Figure  1.  Boundary  Cell  In-  and  Out-Flow  Illustration: 
(a)  and  (a),  also  (b)  and  (E)  Are  Half-Periods  Apart 
(Shaded  region  is  land  area.) 


(Q  The  G ridded  Ocean  Basin  and  Bathymetry.  The  definition  of  a finite-difference  analog  in  place  of  the 
equations  of  oceanic  tidal  motion  requires  a mesh  discretization  of  the  ocean  area  under  investigation.  While  one 
would  like  to  utilize  a grid  system  as  finely  detailed  as  possible  to  achieve  a high  resolution  of  the  tidal  currents, 
one  is,  of  course,  quite  limited  by  the  memory  capacity  of  the  computer  facility  available  and  by  the  known  bathy- 
metric charts  of  the  ocean  basins.  For  the  worldwide  oceans,  essentially  1°  by  1°  averaged  depth  data  were  compiled 
by  Dishon  (1964),  Dishon  and  Heezen  (1968),  and,  independently,  S M.  Smith  et  al.  (1966).  The  latter  collection 
has  been  somewhat  revised  by  Gates  and  Nelson  (1975),  and  both  sets  are  available  in  tape  form. 

In  spite  of  the  specially  collected  depth  data,  all  investigators  prepared  their  own  bathymetric  charts  corre- 
sponding to  the  ocean  grid  system  chosen.  All  data  are  quite  drastically  curtailed  by  artificial  limits  on  both  the 
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shallow  and  deep  ends  of  the  depth  range.  Moreover,  most  computations  were  carried  out  on  rather  coarse  grids  of 
2°  by  2°  to  6°  by  6°.  Obviously,  in  such  an  ocean  geometry,  much  of  the  absolutely  necessary  resolution  of.  for 
example,  narrow  ocean  ridges  and  continental  shelves  is  lost.  In  fact,  all  published  numerically  constructed  tidal 
maps  suffer  from  such  a lack  of  proper  resolution.  One  may  expect  a 1°  by  1°  seamount  or  island  to  have  little 
effect  on  the  surrounding  tidal  currents,  but  a row  of  such  obstacles  over  several  mesh  cells  represents  a hydro- 
dynamical  barrier  separating  the  tidal  motions  on  both  of  its  sides. 

Based  on  the  experience  of  earlier  researchers,  the  author  (Sect.  S.B)  utilizes  a 1°  by  T grid  system  that  is 
"graded"  toward  the  North  and  South  Poles  to  maintain  a more  uniform  mesh  area.  As  was  demonstrated  by 
Zahel  (1970).  such  a grading  is  desirable  to  achieve  a more  uniform  accuracy,  higher  stability,  larger  time  steps, 
and  shorter  computer  times.  The  depth  data  compiled  by  S.  M.  Smith  et  al.  (1966)  were  revised  and  linearly  inter- 
polated to  fit  the  new  grid  system.  At  the  shallow  end  of  the  depth  scale,  the  oceans  were  cut  off  at  the  50-m 
depth  level.  The  new  data  set  retained  a depth  range  from  SO  to  7 700  m,  where  the  upper  limit  remained  un- 
touched. 

The  author's  exploratory  computations  of  the  M , tide  used  the  revised  S.  M.  Smith  et  al.  ( 1966)  depth  data 
and  constructed  a Preliminary  M , Tide  (Schwiderski.  1976)  without  any  empirical  tidal  elevations.  The  overall 
accuracy  attained  was  encouraging.  Over  most  ocean  areas,  satisfactory  agreement  with  numerous  observations  at 
island  and  continental  stations  was  achieved.  Copies  of  this  model  have  been  requested  by  many  research  institu- 
tions for  various  applications.  For  example,  Mr.  Clyde  Goad  at  the  National  Oceanographic  and  Atmospheric  Admin- 
istration's (NOAA)  National  Geodetic  Survey  computed  the  Mytide  effects  on  the  moon's  orbit  and  found,  for 
the  first  time,  good  agreement  with  observations  and  other  theories.  Yet,  a considerable  lack  of  accuracy,  especially 
in  the  phase  constants,  was  encountered  particularly  over  narrow  ocean  ridges,  such  as  the  Aleutian,  Carribbean. 
Marianas,  and  Sunda  Ridges,  paralleling  deep  trenches.  Although  a 1°  by  1 ° grid  system  was  employed,  the  results 
showed  no  satisfactory  improvement  over  computations  made  by  earlier  researchers.  For  example,  the  Pacific  tide 
was  found  sweeping  freely  into  the  Bering  Sea  without  any  separation  effect  expected  from  the  dividing  Aleutian 
Ridge.  This  striking  failure  of  resolution  may  be  seen  in  all  purely  computed  cotidal  maps,  such  as  the  well-scaled 
map  of  Zahel  ( 1 970).  which  also  displays  the  corresponding  observed  values. 

A reinspection  of  the  bathymetric  data  revealed  clearly  that  even  a 1°  by  1°  grid  scheme  falls  far  short  in 
representing  a narrow  ocean  ridge  in  a hydrodynamically  proper  fashion.  The  defect  is  particularly  compounded 
when  the  narrow  ridge  parallels  a deep  trench.  The  reason  for  this  deficiency  is  obviously  the  purely  “hydrostatic" 
character  of  the  averaging  principle  employed  by  S.  M.  Smith  et  al.  (1966)  in  order  to  assign  a depth  value  at  the 
center  of  a mesh  cell  that  is  supposed  to  be  representative  for  the  entire  cell.  For  instance,  if  the  area  of  a mesh 
cell  is  (by  subjective  sight)  more  than  half  land,  then  it  is  called  a “land  cell"  and  the  cell  is  given  (for  the  present 
purpose)  the  depth  value  “zero."  In  the  alternative  case,  the  cell  is  declared  “oceanic"  and  a depth  value  is  assigned 
that  conserves  the  estimated  actual  water  mass.  Because  of  this  “hydrostatic"  principle,  cells  were  found  that  con- 
tained elongated  islands  crossing  even  several  cells,  but  every  cell  was  declared  oceanic.  Moreover,  an  oceanic  trench 
portion  of  the  cell  with  some  7 000-m  true  depth  produced  an  average  depth  of  more  than  3 500  m.  Clearly,  for 
ocean  current  models  the  entire  cell  represents  an  impassable  wall  and  the  depth  value  should  be  "zero"  instead  of 
3 500  m. 

In  order  to  eliminate  the  shortcomings  of  the  bathymetric  tables  compiled  by  S.  M.  Smith  et  al.  (1966)  and  to 
implement  the  necessary  “hydrodynamical"  averaging  principle  recognized  above  (Sect.  S.B).  a renewed  worldwide 
revision  of  the  depth  data  was  undertaken.  More  than  3 000  depth  data  belonging  to  continental,  island,  and  sub- 
merged-ridge cells  were  modified.  In  addition,  the  minimum  and  maximum  cutoff-depth  values  were  set  at  10  and 
7 000  m,  respectively.  The  maximum  limit  affected  only  few  isolated  holes  with  unnecessarily  adverse  consequences 
on  the  stability  (Sect.  5.G)  of  the  tide  model.  Exploratory  computations  with  the  new  ocean-bottom  relief  ini- 
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mediately  reflected  the  improved  resolution.  The  improvement  of  the  results  was  even  more  spectacular  when  empi- 
rical tidal  elevations  (see  (d)  and  Sect.  5.F)  were  hydrodynamically  interpolated. 

Naturally,  it  must  be  kept  in  mind  that  the  decision  ^declare  a boundary  mesh  cell  as  oceanic  and  the  subse- 
quent assignment  of  an  averaged  depth  value  is  (even  with  an  empirically  founded  judgment)  entirely  subjective  and, 
hence,  subject  to  some  error.  While  isolated  inaccuracies  of  this  sort  were  generally  found  ineffective  to  alter  the 
overall  quality  of  the  tide  description,  errors  in  a series  could  not  be  tolerated  in  a satisfactory  tide  model.  Serious 
consequences  of  this  uncertainty  are  greatly  reduced  in  the  present  model,  which  is  supported  by  numerous  tide 
data  observed  at  boundary  stations.  In  fact,  in  the  hydrodynamical  interpolation  procedure  (see  (d)  and  Sect.  5.F), 
it  was  established  necessary  to  allow  some  monitored  periodic  in-  and  out-flow  across  the  mathematically  fixed 
land  side  of  a boundary  cell.  The  unknown,  geographically  true  shoreline  may  be  farther  outside  or  inside  the 
mathematical  cell  and,  thus,  really  require  some  small,  periodic  out-  or  in-flow  (see  Figure  1),  respectively.  Evident- 
ly, the  permissible  in-  or  out-flow  is  enhanced  at  river  estuaries  and  at  narrow  entrances  of  border  seas  (e.g.,  the 
Mediterranean  Sea),  which  are  excluded  from  consideration  as  meshwise  disconnected  from  the  main  oceanic  body. 
The  selected  (Sect.  5.C)  and  incorporated  empirical  tide  data,  however,  do  reflect  the  true  situation  more  closely. 

(g)  The  Finite-Difference  Analog.  On  the  gridded  ocean  basin  (see  (0  and  Sect.  5.A),  the  differential  equations 
of  oceanic  tidal  motions  can  be  converted  into  a discrete,  finite-difference  analog  for  a computerized  solution.  Since 
the  beginning  of  the  computer  age,  numerous  finite-differencing  schemes  with  various  characterisitcs  have  been 
developed  and  tested  by  numerical  analysts  working  in  applied  mathematics,  fluid  mechanics,  oceanography,  meteo- 
rology, etc.  The  application  of  those  techniques  with  or  without  possible  variations  to  the  present  ocean-tide 
problem  depends  essentially  on  the  numerical  tidalist’s  empirically  founded  modeling  skill  and  imagination,  limited 
only  by  the  available  computer  capacity.  Since  the  forcing  function  (see  (c)  and  Sect.  3.D)of  the  basic  tide  model 
is  time-periodic,  the  integration  at  hand  may  be  accomplished  as  a boundary-value  problem  (BVP)  or  as  an  initial- 
value  problem  (1VP). 

In  the  chosen  BVP  approach,  one  removes  the  complex  time-periodic  factor  from  all  three  dependent  vari- 
ables; i.e.,  tidal  elevation  and  east  and  north  velocity  components,  provided  the  basic  tide  model  is  strictly*  linear 
(see  (a),  (b),  and  Sect.  4).  In  the  remaining  complex  (time-independent)  differential  equations  for  the  complex 
amplitudes  of  the  independent  variables,  all  spacial  derivatives  may  then  be  replaced  by  central  finite  differences  cor- 
responding to  the  grid  system.  As  was  first  pointed  out  by  Richardson  (1922),  the  differencing  scheme  may  be 
chosen  “regular,”  “semistaggered,”  or  (fully)  “staggered;”  that  is,  the  three  dependent  variables  may  be  tabulated 
(computed)  at  collocated,  semidislocated,  or  fully  dislocated  points,  respectively,  as  shown  in  Figure  2. 

The  resulting  three  difference  equations  may  be  combined  into  two  equations  for  the  complex  velocity  ampli- 
tudes by  eliminating  the  complex  tidal  amplitude,  which  is  directly  computable  from  the  velocity  via  the  conti- 
nuity equation.  With  sufficient  velocity  boundary  conditions  (no-flow  and,  if  necessary,  free-slip  or  no-slip;  see  (d)), 
one  winds  up  with  an  implicit,  complex  system  of  linear  equations  that  can  be  solved  by  modem  computerized 
methods.  If  higher-order,  eddy-dissipation  terms  are  included  in  the  basic  tide  model,  the  linear  system  to  be  solved 
becomes  quite  involved.  Such  a system  has  not  yet  been  tried.  However,  if  eddy  dissipation  is  completely  neglected 
(see  (a)  and  Sect.  3.C),  the  linear  system  to  be  solved  is  considerably  simpler  and  has  been  successfully  inverted  by 
Accad  and  Pekeris  (1963)  and  Pekeris  and  Accad  (1969). 

The  same  three  complex-difference  equations  derived  above  for  the  three  dependent  variables  can  also  be  com- 
bined into  a single  difference  equation  for  the  complex  tidal  amplitude,  provided  all  eddy -dissipation  terms  are 


•If  nonlinear  properties  (e.g.,  square  law  of  bottom  friction)  are  modeled,  infinite  Fourier  scries  could,  in  principle,  be  used,  but 
those  lead  to  an  involved  infinite  system  of  equations. 
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Figure  2.  Finite-Difference  Schemes  in  Space:  (a)  Regular, 
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neglected.  This  alternative  approach  to  the  linear  BVP  is  obviously  advantageous  whenever  empirical  tidal-boundary 
values  are  introduced  in  place  of  the  purely  mathematical  velocity-boundary  conditions  imposed  in  the  first  ap- 
proach. Of  course,  the  obtained  simpler-difference  analog  must  still  be  integrated  by  standard  computer  programs. 
This  oversimplified  model  attracted  the  attention  of  several  investigators,  such  as  Hansen  (1948,  1949),  Bogdanov 
and  Magarek  (1967, 1969),  and  Hendershott  (1972). 

Since  the  basic  tide  model  is  time -dependent,  it  may  be  integrated  as  an  IVP  starting  from  some  specified 
initial  state.  In  addition  to  the  advantages  already  pointed  out  under  (d)  above,  this  most  natural  approach  leaves 
unrestricted  room  for  full-scale,  realistic  modeling  of  oceanic  tidal  currents,  including  linear  or  nonlinear  eddy 
dissipation  (a)  and  bottom  friction  (b)  as  well  as  mathematical  boundary  data  together  with  hydrodynamically 
interpolated  empirical  tidal  data  (d).  In  fact,  in  Section  6.A  it  will  be  argued  that  it  is  the  discrete,  time-dependent, 
finite-difference  analog  that  reflects  the  physically  realistic  properties  of  oceanic  tidal  currents  much  more  perfectly 
and  ponderably  than  the  continuous  differential  model.  The  differential  equations  can  be  used  advantageously  to 
define  various  finite-difference  formulations,  but  it  is  neither  possible  nor  desirable  to  seek  convergence  of  the 
discrete  solution  to  the  continuous  integral  as  is  the  objective  in  laminar-flow  problems. 

The  finite-difference  scheme  of  the  BVP  discussed  above  remains  the  same  for  the  IVB  in  space  variables.  In 
the  differential  equations  of  oceanic  tidal  motion,  derivatives  in  space  variables  may  be  replaced  by  central  finite 
differences  employing  a Richardson  regular,  semistaggered,  or  (fully)  staggered  scheme.  The  resulting  equations 
conta  h first-order  derivatives  in  time,  which  can  be  inverted  by  integration  over  one  or  two  time  steps  along  the 
degenerated  characteristic*  of  the  system,  i.e.,  parallel  to  the  t-axis.  The  remaining  (not  explicitly  integrable)  in- 
tegrals are  then  replaced  by  some  average  integration  rule  (e.g.,  one-point  tangential,  two-point  trapezoidal,  three- 
point  Simpson).  Depending  on  the  chosen  time-integration  formula,  the  obtained  finite-difference  analog  must  be 
solved  by  a computer  applying  either  an  explicit  or  implicit  time-stepping  procedure  that  starts  from  an  initial  state 
and  incorporates  mathematical  (no-flow,  free-slip,  or  no-slip)  and  empirical  boundary  data  (see  (d)). 


•It  may  be  suggested  here  that  an  integration  along  the  ruled  lines  of  the  characteristic  cone  of  the  system  might  be  preferable;  but 
this  alternative,  or  some  mean  of  both  possibilities,  has  never  been  tried  in  the  modeling  of  ocean  currents. 
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As  was  recognized  by  Richardson  (1922),  the  time  integration  may  utilize  a “regular,”  a “step-over”  (“leap- 
frog”), or  an  “intermediary”  (“alternating”)  time-stepping  scheme.  In  the  “regular”  scheme,  all  three  (real) 
dependent  variables  (east  and  north  velocity  components  and  tidal  height)  are  computed  and  tabulated  at  the 
same  time  points  (Figure  3a).  In  Richardson’s  “step-over”  (now  often  called  “leap-frog”)  method,  both  velocity 
components  are  co-timed,  but  tidal  elevations  are  computed  and  tabulated  at  intermediate  points  (Figure  3b). 
In  the  “intermediary”  (or  alternating)  procedure,  all  three  dependent  variables  are  computed  and  tabulated  coin- 
cidently,  but  both  velocity  components  use  computationally  auxiliary  intermediate  time  points  in  an  alternating 
fashion  (Figure  3c). 
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Figure  3.  Finite-Difference  Schemes  in  Time:  (a)  Regular,  (b)  Step-Over  (Leap-Frog), 
(c)  Intermediary  (Alternating);  (U.  V.  J)  Computed  and  Tabulated, 

(V,  V)  Computed  Not  Tabulated 


Although  the  step-over  scheme  is  widely  used  in  general  ocean-current  and  atmospheric  circulation  studies, 
the  regular  and  intermediary  schemes  seem  to  have  been  preferred  in  tidal  models.  The  intermediary  method  has 
been  applied  in  a rather  involved  explicit-implicit  form  by  Marchuk  et  al.  (1969),  Gordeyev  et  al.  (1973),  and 
Marchuk  et  al.  (1973)  to  tidal  computations  in  border  seas  and  in  a global  ocean  model.  In  the  now  well-known 
hydrodynamical  numerical  method  of  Hansen  (1966),  the  staggered  scheme  in  space  is  combined  with  the  regular 
scheme  in  time  in  an  explicit  and  convenient  fashion.  The  method  was  successfully  applied  by  Zahel  (1970, 1973, 
1975)  and  Estes  (1975,  1977)  to  the  tides  in  global  oceans  on  4°,  3°,  2°,  and  1°  square-grid  systems. 

The  principles  of  the  Hansen  technique  have  been  selected  as  most  appropriate  for  the  present  investigation 
because  of  their  simplicity  and  versatility.  As  will  be  shown  in  Section  5.E?  the  staggered  scheme  in  space  permits  a 
very  easy  work-in  of  no-flow  and  free-slip  (or  no-slip)  boundary  data.  The  integration  rule  over  one  time  step  used 
by  Hansen  will  be  improved  in  Section  5.D  by  using  mixed  averages  without  losing  the  important  explicit  property 
of  the  finite-difference  scheme.  This  modification  resulted  in  an  enhanced  stability  of  the  system.  Moreover,  it 
facilitated  the  implicit  hydrodynamical  interpolation  of  empirical  tide  data. 
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3.  THE  CONTINUOUS  OCEAN  TIDAL  EQUATIONS  (COTES) 


A.  THE  NAVIER-STOKES  EQUATIONS  OF  AVERAGED  TURBULENT  FLOW 

Because  of  the  enormous  dimensions  of  the  world  oceans,  their  hydrodynamical  tidal  motion,  generated  by 
the  attraction  of  the  moon  and/or  sun,  must  be  considered  entirely  turbulent.  Accordingly,  any  comprehensive 
modeling  of  oceanic  tidal  currents  should  begin  with  the  complete  Navier-Stokes  equations  of  averaged  turbulent 
motions  of  a viscous,  incompressible  fluid  including  all  unknown  Reynolds  stresses  (see,  e.g.,  Schlichting,  1968), 
which  must  provide  the  major  friction  stress  to  keep  the  intrinsic  supercritical  instability  of  the  flow  under  control. 
In  rotating  spherical  polar  coordinates  (X,  0.  r).  these  equations  may  be  written  in  the  form  (see,  e.g.,  Whitaker, 
1968) 
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In  these  equations,  all  subscripts  denote  indicated  partial  derivatives,  while  all  superscripts  denote  indicated  com- 
ponents of  the  turbulent  dissipation  vector  f or  the  Reynolds  stress  tensor  r.  Furthermore,  the  following  notations 
are  used: 

t = universal  time  in  sec 


= east  longitude 

= _ _ s = north  latitude  ( 0 = colatitude) 
= polar  radius  in  m 
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(//,  v,  n)  = east,  north,  and  radial  velocities,  respectively,  in  m/sec  (averaged) 

/>  = pressure  (averaged) 

</  = total  tide-generating  potential 

/ = dissipation  vector 

r = Reynolds  stress  tensor  (to  be  specified) 

S2  = 0.72722  X 10'4  sec'1  = earth  angular  velocity 

p 103  kg/m3  = density  of  sea  water 

(»'  = 9.8 1 m/sec‘  = gravity  acceleration 

The  equations  of  turbulent  "mean"  motion  are  obtained  by  a formal  time-averaging  procedure  applied  to  the 
Navier-Stokes  equations  of  viscous  laminar  flow  that  leads  to  the  concepts  of  “averaged"  velocities  and  pressure, 
as  well  as  to  the  unknown  Reynolds  stress  tensor,  r.  containing  the  filtered  out,  fluctuating  velocity  residuals  in 
quadratic  form  (see,  e.g.,  Schlichting,  1968).  A discussion  of  the  physical  meaning  of  these  and  following  tur- 
bulence notions  may  be  postponed  to  Section  6. A.  The  momentum  equations  (Eq's.  1)  maintain  the  momentum 
balance  between  the  familiar  kinetic  (inertial  acceleration)  forces  on  the  left  side  and  the  forces  of  potential  (tide, 
gravity,  and  pressure),  Coriolis,  centrifugal  acceleration,  and  dissipation  of  the  right.  The  continuity  equation 
(Eq.  2)  expresses  the  condition  of  conservation  of  mass.  The  Reynolds  stress  tensor,  t,  and  the  total  tide-generating 
potential.  </.  will  be  specified  in  Sections  3.C  and  3.D. 


B BOTTOM  AND  SURFACE  BOUNDARY  CONDITIONS 
The  global  oceans  at  hydrostatic  rest  may  be  described  by 

(а)  r = R = 0.637  X I07m  = spherical  (geoidal  rest)  sea  surface  (all  geoidal  undulations  are  neglected  without 
any  significant  loss  of  accuracy). 

(0)  r = R //(X.  0)  = sea-bottom  relief,  where 

(y)  //  = //(X.  0)  = realistic  ocean  depth  in  m (//(X,  0)  = 0 for  land  points,  see  Sect.  5.B).  and 

(б ) p = P (ip/.  = hydrostatic  sea-pressure  distribution  with  constant  (arbitrary)  sea-surface  pressure  P.  where 
(e)  z=r  R=  new  depth  variable,  so  that  z - 0 denotes  r =R  (see  Figure  4). 

Due  to  the  time-dependent,  tide-generating  potential,  q.  acting  on  the  ocean  and  solid  earth,  the  hydrostatic 
conditions  (a)  (e)  are  altered  to  the  following  hydrodynamical  boundary  conditions  of  the  sea  surface  and  bottom 

(see  Figure  4): 

(a)  zs  = fs  = fs(X.  0.  t ) = total  sea  surface  tidal  elevation  over  geoid  z = 0; 

(b)  z*  = f*(X.  <)>.  ()  //(X.  0)  = sea-bottom  tidal  relief,  where 

(c)  = f^’lX.  0.  t)  = total  bottom  tidal  elevation  over  rest  relief  z = H(\.  0); 

(d)  f = f(X.  0.  M = Js  f7’  = ocean  tidal  elevation  (measured  by  bottom  tidal  pressure  gauges)  to  be  modeled: 

(e)  if  = r<X.  0. 1)  = earth  tidal  elevation  to  be  specified  (Sect.  3.D.  Fq’s.  9); 

(0  = ff,’(X.  0.  r)  = f = earth  dip-response  to  oceanic  tidal  load.  J (Eq.  1 0): 

(g)  ps  =P  = constant  surface  pressure  (no  atmospheric  pressure  considered): 

(h)  w*  = f!  = surface  radial  velocity; 

(1)  T5  = o = no  surface  stress  (free-slip.  no  wind  force  considered): 

(j)  (uh.  vh.  »h)  • V(//  + -)  ^ = no  flow  across  ocean  bottom  (V=  gradient  vector): 

(k)  rh  = pB(uh.  vh)  = bottom  stress  vector  specified  by  invoking  the  linear  law  of  bottom  friction  (Sect.  2(b)) 
with  the  coefficient 

H = bp  cos  0 = hi.  ’/ucos  0 ( 4a ) 
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Z (r) 
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Z*  * f*  * TIDAL  SURFACE 
Z * 0 * GEOIDAL  SURFACE 

Zb  = Cb-H* BOTTOM  TIDAL  RELIEF 
Z = -H  = BOTTOM  REST  RELIEF 


Figure  4.  Earth-Ocean  Tidal  Interaction:  f = Ocean  Tide,  ff  = Earth  Tide,  = Surface  Tide, 
= Bottom  Tide,  = if  - p = Earth  Dip-Response  to  Ocean  Tide, 
and  //  = Ocean  Depth 


depending  on  the  cell  bottom  area  L2f jcos  0(Sect.  6. A),  where 

(l)  L = chosen  equatorial  mesh  size  (Sect.  5. A), 

(m)  = mesh  “grading"  parameter  (Sect.  5. A),  and 

(n)  b = bl  2 = uniform  bottom  friction  parameter,  which  must  be  determined  by  trial-and-error  computations 
for  best  results.  The  Mj  tide  (Sect.  6.B  and  Part  II)  was  computed  with 

b = 0.01  m/sec,  (4b) 

which  is  considerably  below  the  value  used  by  Pekeris  and  Accad  (1969).  Some  implicit  local  adjustments  of  b will 
be  allowed  in  Section  S.F  to  better  accommodate  observed  tidal  data. 

(o)  Lateral  boundary  values  will  be  specified  in  Section  S.E. 


C.  REYNOLDS  STRESSES  AND  EDDY  DISSIPATION 


The  formally  derived,  unknown  Reynolds  stress  tensor,  r (Sect.  3. A),  may  be  specified  by  the  original  Bous- 
sinesq  (1877)  substitution  (Sect's.  2(a)  and  6. A);  that  is,  by  the  laminar,  viscous  stress  tensor  (see,  e.g.,  Schlichting, 
1968;  Whitaker.  1968): 
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In  this  substitution,  the  ordinary  kinetic  molecular  viscosity  is  replaced  by  the  so-called  “eddy  viscosity"  (momen- 
tum austausch  exchange,  mixing  coefficient),  which  remains  to  be  modeled  to  represent  the  true  turbulent  flow 
characteristics  at  hand  as  closely  as  possible. 


In  order  to  achieve  a greater  modeling  flexibility,  it  is  customary  to  divide  the  eddy  viscosity  into  a “vertical" 
eddy  viscosity  associated  with  vertical  shear  and  into  a "horizontal"  (lateral)  eddy  viscosity*  associated  with  hori- 
zontal shear.  In  single-layer  ocean  models  like  the  present  one  (Sect.  3.E).  the  separate  treatment  of  the  vertical 
eddy  viscosity  needs  no  explicit  specification,  because  (Eq’s.  28)  the  viscosity  becomes  an  important  part  of  the 
bottom-friction  coefficient.  B.  defined  by  Equations  4a  and  4b. 


Based  on  the  author's  extensive  computer  experiments  and  on  the  physical  arguments  elaborated  in  Sections 
2(a).  5.G.  and  6. A the  (horizontal)  eddy  viscosity.  A . may  now  he  specified  by 

A =^/.//(X.0KI  +/jcos  <P).  (6a) 

In  this  definition.  /.  and  p denote  the  mesh  size  and  grading  parameters  introduced  in  Section  3.B(I)  and  (m).  Ac- 
cordingly. the  novel  eddy  viscosity  depends  on  the  mean  lateral  cross-section  area  //  X /.( I + pcos  0)/2  of  the  flow 
cell  of  depth  H and  average  northeast  mesh  size  /.(l  + pcos  <p)/2.  The  remaining  reduced  eddy-dissipation  coef- 
ficient, a (in  sec),  must  be  subjected  to  trial-and-error  computations  in  order  to  achieve  best  results  uniformly**  over 


•Since  east  and  north  dimensions  are  usually  nearly  equal,  no  need  lor  two  lateral  eddy  viscosities  has  ever  been  encountered. 
••Different  a-values  for  the  Pacific.  North  and  South  Atlantic,  and  Indian  Oceans  were  also  tested  without  significant  effects. 
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all  oceans.  It  may  be  recalled  from  Section  2(a)  that  some  dependence  of  A on  the  mesh  size  L was  noticed  by 
several  earlier  researchers.  Clearly,  the  strong  variations  of  L and  the  depth  H explain  the  huge  magnitude  differences 
of  the  “constant”  eddy  viscosity  used  by  analysts  in  different  problems.  In  the  present  case  with  H varying  from 
10  to  7 000  m (Eq’s.  50,  51,  71,  103,  and  123),  the  eddy  viscosity  varied  between 

1.3  • 103—  <A  < 1.3  • 106— , (6b) 

sec  sec 

which  fits  the  customary  range  very  well.  It  is  perhaps  fortuitous  yet  interesting  to  note  that  for  a cubic  cell  of  about 
1-in.  (2.54-cm.)  mesh  size  the  eddy  viscosity,  A (Eq’s.  6a,  b),  reduces  to  the  value  of  the  molecular  viscosity  of 
water. 

D.  THE  TOTAL  TIDE-GENERATING  POTENTIAL 


The  total  tide-producing  potential,  q (Sect.  3. A),  may  be  expressed  in  the  form 

<7=C(W),  (7) 

where  Gtj  is  the  “primary”  astronomical  potential  directly  proportional  to  Newton’s  equilibrium  tide  rj  (see  Sect.  2). 
The  remaining  “secondary”  potential  Gtj’  can  be  expanded  into  its  three  major  parts 

Gtj’  = G(j?°  + tj£  - rjco),  (8) 

which  reveal  their  corresponding  origins  (Figure  4,  Sect.  3.B(d),  (e),  (0.  respectively); 

Grp  = gravity  potential  perturbation  due  to  the  ocean  tide  f 
Gif  = gravity  potential  perturbation  due  to  the  earth  tide 

Gtf°  = gravity  potential  perturbation  due  to  the  earth  dip-response  >f°  caused  by  the  ocean  load  (tide)  f 


Attention  to  the  significance  of  the  earth-ocean  tidal  interactions  manifested  in  the  five  quantities  rf , 

if°,  rf° , and  tj°  was  called  for  by  Farrell  (1972  a,  b;  1973).  As  was  argued  in  Section  2(c),  it  is  sufficient  in  ocean- 
tide  models  to  use  the  following  approximate  relations: 


Sf  *=  0.61tj,  rf  * 0.30tj, 

(9a) 

f - rf  0.3 It?, 

(9b) 

-rjeo  +tj°  *=0.10f. 

(10) 

The  first  relations  (Eq’s.  9)  used  by  Hendershott  (1972,  1975),  Zahel  (1975),  Estes  (1977).  and  others,  can 
be  justified  by  second-order  Love  number  approximations.  The  second  relation  (Eq.  10)  has  been  suggested  by 
Pekeris  (see  Sect.  2(c)),  who  also  recommended  the  factor  0.10  after  evaluating  the  Green’s  function  representation 
of  the  three  oceanic  tidal  load  effects  (f*0,  rf°,  rf)  derived  by  Farrell.  Apparently,  the  suggestion  by  Pekeris 
(Eq.  10)  is  physically  just  as  plausible  as  the  accepted  approximation  (Eq’s.  9).  In  Equations  9a  and  9b,  it  is  assumed 
that  the  solid-earth  tide,  $*,  and  its  subsequent  gravity  perturbation,  Grf , are  essentially  instantaneous  responses  to 
the  tide-generating  potential,  q.  Analogously,  in  Equation  10  it  is  assumed  that  the  solid-earth  dip,  if° , and  the 
gravity  perturbations,  r\eo  and  rf , are  almost  instantaneous  responses  to  the  ocean’s  tidal  load,  f.  It  may  be  men- 
tioned that  the  author  conducted  extensive  computer  experiments  using  the  factors  0.00,  0.08,  and  0.12  instead 
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of  0.10  in  Equation  10.  The  last  two  factors  produced  no  noteworthy  alterations  of  the  results.  The  first  factor 
(0.00)  obviously  deletes  all  oceanic  tidal-load  effects  to  which  critically  large  significance  has  been  attached  by 
Hendershott  (1972,  1975).  The  author’s  computations  supported  the  marginal  effects  of  oceanic  tidal  loading 
found  by  Estes  (1977). 

Following  Thomson  (Lord  Kelvin,  1868),  G.  H.  Darwin  (1883),  Doodson  (1921),  and  Cartwright  and  Tayler 
(1970),  the  primary  astronomical  tide-generating  potential,  6 ),  or,  equivalently,  the  equilibrium  tide,  tj,  may  be 
expanded  into  a series  (see,  e.g.,  Dietrich,  1963)  of  “harmonic  components”  (constituents),  t\v,  with  a nonharmonic 
frequency  spectrum 

^ = 2 n (X,  0,  f).  (11) 

*>=0 

In  the  final  analysis  of  ocean  tides,  the  following  three  major  species  will  have  to  be  considered: 

(a)  Semidiurnal  Equilibrium  Tides 

v = 2:  rj2  =Kcos2  <pcos(ot  + 2X  + x)  (12) 

(b)  Diurnal  Equilibrium  Tides 

v = 1:  ijj  =Ksin  20cos(or  + X + x)  (13) 

(c)  Long-Period  Equilibrium  Tides 

v = 0:  r?Q  =^1  - 3cos2  0)cos(ar  + x)  (14) 

In  Equations  12-14,  the  constants  ( K , a,  x)  denote 
K = amplitude  of  equilibrium  tide  (in  m) 
a = frequency  of  equilibrium  tide  (in  sec'1) 

X = astronomical  argument  of  equilibrium  tide  (in  rad) 

In  Table  1,  these  constants  are  listed  for  the  major  tidal  modes  with  amplitudes  larger  than  4 percent  of  the 
leading  semidiurnal  moon  (M2)  tide.  The  daily  astronomical  argument,  x.  can  be  neglected  in  the  following  construc- 
tion of  the  oceanic  tidal  modes 


f = £(X,  0 ) cos (at  + x - 6(X.  0))  (15a) 

corresponding  to  the  considered  mode  of  the  equilibrium  tide,  77  =V„.  According  to  Equation  15a,  only  the  “har- 
monic constants” 

(i  = 0)  * tidal  amplitude  (in  m)  1 


and 


(15b) 


6 = 5 (X,  0)  = Greenwich  phase  (in  rad) 


need  to  be  found,  provided  the  tide  model  is  linear  or  almost  linear  (Sect.  2(g)). 
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E.  SIMPLIFYING  ASSUMPTIONS 


In  addition  to  the  simplifying  assumptions  made  concerning  bottom  friction  (Sect.  3.B(k)-(n)),  eddy  dissipa- 
tion (Sect.  3.C),  and  total  tide-generating  potential  (Sect.  3.D),  the  following  simplifications  may  be  invoked: 

(a)  Hydrostatic  pressure  assumption  in  Equations  1 ; i.e., 

(а)  Neglect  all  quadratic  inertial  accelerations  (Stokes  slow-motion  assumption  consistent  with  linear 
eddy  dissipation  and  bottom  friction  (see,  e.g.,  Schlichting,  1968)) 

(0)  Neglect  all  centrifugal  accelerations 

(>)  Neglect  vertical  Coriolis  force 

(б)  Neglect  vertical  dissipation 

(e)  Neglect  vertical  motion 

(b)  Single-layer  ocean  assumption  in  Equations  1 and  2;  i.e., 

(a)  r = R + z * R,  but  dr  = 9z 

(0)  f(X,0,r)<</f(X,0) 
f*(X,  0,  t)<H(\,<t>) 

(y)  u(X,  <t>,  z,  r)  * u(X,  0,  r) 
v(X,  0,  z,  f)  ^v(X,  0,  t) 

It  may  be  pointed  out  that  all  assumptions  of  (a)  and  (b)  are  well  justified  over  most  flow  areas.  In  particular, 
the  strong  assumptions  ( b , y)  are  realistic  because  the  tide-generating  potential  is  a body  force.  Moreover,  the  motion 
is  fully  turbulent  and,  hence,  the  averaged  velocity  profile  exhibits  only  a very  thin  boundary  layer  and  laminar 
sublayer  (see,  e.g.,  Schlichting,  1968).  It  may  be  mentioned  that  for  the  same  reason  the  condition  of  free-slip  at 
a boundary  wall  (Sect’s.  2(d)  and  5.E)  seems  more  appropriate  than  the  no-slip  condition  used  in  laminar-flow 
situations.  The  assumption  ( b , y)  is,  of  course,  much  less  realistic  in  general  ocean  currents,  which  are  driven  by 
surface  pressure  and  wind  forces  and/or  by  density  variations  confined  to  the  upper  ocean  layers. 

By  applying  assumption  (a)  to  Equation  lc,  one  finds  the  hydrostatic  pressure 


or,  with  (Sect.  3.B  (a)  - (g)) 


P = P , 


rs  = * + = i + r - r°. 

P = Gp(!  + $e  - se°  - z)  + P . 

With  Equations  7-10  and  16,  one  has 

q - = G [17  - (ie  - r,e)  - r + (ff°  - r,eo  + 170)]  - />/p, 


q G (orr?  - 0f)  - P/p. 


Where, 


a = 0.69, 


0 = 0.90. 


In  the  following  equations,  it  is  convenient  to  introduce  the  notations 
"a  = "JH-  "ax  - HJ«- 

"♦  - HJH-  - HJH- 


H = JT.  - V sin  0/(1  + /i  cos  0) 

0 0 

so  that  (Eq.  6a) 

/4  = A ffx  , A B A H . <21> 

X X 0 0 

Using  the  simplifications  (a)  and  (b)  and  Equations  17-21  in  Equations  1 and  2,  one  arrives  at  the  reduced 
equations  of  motion: 

U,  = R (0TJ  ' + 2«vsin0  ♦ f\  (22a) 


V'  = R ^a,?  " - 2fl  m sin  0 + /*, 
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and 


mx  + (v  cos  0)^  + R cos  0 w2  = 0. 
where  (Eq’s.  3a,  b;  5a-f;  and  2 1 ) 

Ri  £>0  Kx  + 2">x  - “1  + A Uzz 

+ Ko  + (^0  ~ ,an^“0  + V ‘an  0] 


+ -/F^?0  l(W0  “ 2tan0)vx  - 2##xvtan0], 


/»=  K,  + //.v4  - v]  + Av 

J R 1 cosJ  0 1 xx  xx  1 z 


+ jp  K* + (2//0  _ ,an  *)v<»] 


+ 'RT^0  I(2mx  + *x“),a"*  + "x“J- 


F.  DERIVATION  OF  CONTINUOUS  OCEAN  TIDAL  EQUATIONS 

Under  the  single-layer  ocean  assumptions  made  in  Section  3.E(b),  the  reduced  Equations  22  and  23  may  be 
integrated  over  the  instantaneous  ocean  depth  (zJ  - zb  = H + f ) while  observing  the  surface-  and  bottom-boundary 
data  specified  in  Section  3.B.  For  that  purpose,  it  is  useful  to  introduce  the  “integrated”  velocity  components: 

fZS 

l/(X,  0,  i)  = I udz  *uH  = (integrated)  east  velocity  (25a) 

Jzb 

and  ^ 

V(X,0,  t)=  I '/dz  % = (integrated)  north  velocity,  (25b) 

J2b 

where  the  term  “integrated"  may  be  omitted  when  no  confusion  appears  possible. 

Using  the  simplifications  and  notations  of  Section  3.E,  one  finds  the  helpful  approximations  with  the  corre- 
sponding replacements  (m  **  v,  U **  V,  X ♦*  0): 


,z' 

J u(dz  mUf  — u*z*  + ubzbt  » U(, 

*2b 

fzs 

I k uxdz  = Ux  - u‘<  + “*zx  * Ux  - W 


and 


/„  «XX*  * "xX  "x"  + ("l  ~ ~ “A'  (26C) 

The  bottom-boundary  conditions  imposed  in  Section  3.B(j)  and  (k)  assume  the  following  approximate  forms: 

(27) 

(28a) 

(28b) 

With  Equations  25-28  and  the  surface-boundary  conditions  in  Section  3.B(h)  and  (i),  the  proposed  integration 
of  Equations  22  and  23  is  easily  carried  out  and  yields  the  following  “continuous  ocean  tidal  equations”  (COTEs 
with  9 = rr/2  - <p  = colatitude): 


f*  = — 
s t R 


^ub  + v/ 


cos  <p 


i T*x  * a(  ^ + «A 

P T H y?2cos20  R 2 7 


and 


1 fc*  B ( HA  "A  rb) 
_T  --  + r2  , ) ■ 


+ £ |"xx  + "x^x 


/?2  l sin2  9 
A 


+ l/„  + (cot  9 - + We)f/9 


/?2  sin  0 


[HKHeV  - (2  cot  9 + //„)KX]  + 2«Kcos0, 


1 ‘ -f  « - - K(f  ♦ F "*) 

♦ £ pi-  * «w  * (»'«  - «» * 


R2  sin  9 


[2  cot  9 - HxUe  - (cot  9 - He)HxU]  - 2Ot/cos0, 


and 


(29a) 


(29b) 


R sin  9 — (K  sin  9)e  = 0, 


(30) 
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where  a = 0.69,  0 = 0.10,  and 


x + \ + Hi  _ __ 

H sin2  e — ^ + Hee  ~ He(H*  " He)  + (H$  * H°)cotd’ 

Hxx  + 1 

H»=  + Hee  + He(cot0  - He  + 2 He), 

He=  Hg/H,  Heg  = Hee/H. 

He-  He  + ji  cos  0/(1  + n sin  0), 

A = LH(\<9)( l+jtsin0),  B = 6/isin0. 


(31a) 

(31b) 


(32) 


(33) 


It  may  be  mentioned  that  for  a = 0 = 1 and  A = B = 0 the  ocean-tide  equations  (Eq’s.  29  and  30)  reduce 
to  the  considerably  simpler  classical  Laplace  tidal  equations.  Evidently,  the  complete  COTEs  require  second  deriva- 
tives (Hxx,  Hee ) of  the  bottom  topography,  which  can  be  assumed  to  exist  without  placing  major  restrictions  on 
the  realistic  features  of  the  bathymetry.  It  is  interesting  to  note  that  those  second  derivatives  in  H x and  H6  act 
as  modifying  bottom-friction  terms.  A ridge-like  ocean  floor  (say  ,HXX  > 0)  always  adds  to  the  bottom-friction  terms 
(U.  V)  B/H , while  a valley  (Hxx  < 0)  always  diminishes  bottom  friction. 
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4.  SECOND-ORDER  ARCTIC  TIDES 

As  was  mentioned  in  Section  2(e),  the  COTEs  (Eq’s.  29  and  30)  become  singular  at  the  North  Pole.  For  the 
intended  numerical  procedure  it  is  therefore  advantageous  to  seek  an  approximate  analytic  solution  around  this 
singularity  that  can  be  matched  together  with  the  numerical  solution  at  some  appropriate  colatitude.  In  fact,  a 
unique  “second-order  arctic  tide”  solution  can  be  determined  for  all  three  species  of  tide-generating  potentials 
listed  as  Equations  12,  13,  and  14,  provided  the  ocean  depth  around  the  North  Pole  is  assumed  constant. 

The  COTEs  (Eq’s.  29  and  30)  for  constant  depth,  H = HQ,  constant  eddy  viscosity,  A = AQ,  and  constant 
bottom-friction  coefficient,  B = BQ,  assume  the  simpler  form 

+ pi  A°2  a WK.-U+  sin  0 (sin  6 Ue)e  - 2 cos  0 FJ,  (34a) 

/\  sin  u 

,,  GH0  B0 

Vt  = ~r~  (fit  ~ «»?)«,  - 2ni/cos0  -~V 

+ D1  A°2  0 — K + sin  0 (sin  6 Ve)e  +2  cos  0 t/J , (34b) 

i\  sin  u 

and 

f'  + A1HT0  [U*  ~ (Ksin0)o)  = °-  05) 

The  forcing  equilibrium  tides,  rj  (Eq’s  12, 13,  and  14),  may  be  written  in  the  unified  complex  form 

17  = rT(X,0)e'a',  (36) 

where 

K sin*  0 e2iK  , (p  = 2)  (37a) 

r?  =r?p(M)=  A sin  20  e,x  , = 1 ) (37b) 

j(l  - 3sinJ  0),(«'  = O).  (37c) 

With  the  substitution 

(v,(.  v.i 0 = (n,l, 


one  arrives  at  the  following  three  complex  differential  equations  in  (X,0): 


aU  = -b'H° a (P?  - arl)\  + 2S2K  cos  0 + i~  U 
sin  u **o 

- [{7X.  - U + sin  0 (sin  9 De)e  - 2/  cos  9 Px| , 

R sin2  0 

a?  = 9Ss  (aq  - fif).  + 2SlUcos9  + / -£-  V 

R n° 

_ — — — [Fx,  - V + sin  0 (sin  9 Ve)e  - 2i  cos  0 £/.  ] , 
R 2 sin2  0 

°f=  + (Psin0)«]‘ 


(39a) 


(39b) 


(40) 


The  reduced  Equations  39  and  40  may  be  solved  for  regular  solutions  by  power  series  in  sin  9 (essentially 
polar  distance)  of  the  form : 


For  v = 2 and  v = 0, 

F = fo(X)  + f i (X)  sin  9 + f2(X)sinJ0  + ..., 

U = t/0(X)  + U i (X)  sin  0 + f/2(X)  sin20  + ..., 

and 

V = cos  0 IK0(X)  + F ! (X)  sin  0 + F2(X)sin20  + ] ; 

and  for  v = 1 , 

r = cos 0 (f0(X)  + f i (X) sin 0 + ?a(X) sin20  + ...], 

U = cos  0 [f/o(X)  + Ui  (X)  sin  0 + t/2(X)  sin2  0 + . . .] , 
and 

F = F0(X)  + F,(X)  sin  0 + F2(X)sin20  + ...  . 


(41) 


(42) 


By  truncating  the  series  expansions  (Eq’s.  41  and  42)  after  the  second  order  in  sin  0 and  substituting  these 
truncations  into  Equations  39  and  40  up  to  the  same  quadratic  power,  one  arrives  after  some  lengthy  but  simple 
algebra  at  the  following  unique  “second-order  arctic  tidal  approximations.” 


For  semidiurnal  tides  (v  = 2), 


r>  = K sin20  e'(2x  + a,\ 
f = 6^sin30e'(2A  + o,), 

U = -2KoRsu\6  eK2x* 
and 

V = -iKoR  sin  6 e,(2x  + 0,) , 

where 

K = uGHoK/[60GHo  + o/?2(2S2  - o)  + ra(6/40  + B0R2/H0)]\ 

and  for  diurnal  tides  (v  = 1), 

tj  = K sin  26  e**A  + 0,) , 
f = 3K,  sin20e'(A  + a'\ 

U = (£j  + 6A:3sin20)cos0e'(A  + o')  , 
and 

K = i [tf,  + 2(^3  + K,0/?)sin20]e'(A  + o,), 

where 


(43) 


(43a) 


(44) 


! 


_ ^31^33(^22  —^12) 

3 R22  (^13^31  - Ri\K3i)  - Ki2(K3iKi3  - K21K32) 

^2  = P 1^33  - ^3(^13  - ^11^3j/^3l)]. 

* 12 

Ky  = —K  3^(32/^131  ; 


(44a) 
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and 


/Cn  = (60GHO  - a2R2)  + io(6A0  + B0R2/H0 ), 
Kl2  = -SIR, 


Kl3  = R(6Sl  - a)  + ^(1240  + B0R2/H0), 

*2,  = 60GHO  + 4<a40  , 

*22  = K(o  - 2SI)  ~^(2A0  + B0R2/H0), 

K23  = \6iA0/R , 
tf31  = 2o£lR2 , 

ATjj  = R(2Sl  - 3a)  + 3^(1240  + B0R2/H0), 

K33  = 2aGH0K. 

For  long-period  tides  ( v = 0), 

tj  = 4"  0 -3  sin^e'0*  , 
f = 2K,  (2  - 3 sin30)  e,of , 

U = K2  sin  6 e'°l , 
and 

V = iK  yoR  sin  20  eial  , 

where 

£3  = 4oSlR3l[oR 3 - i (14 0 + fio/?J/^o)l, 

*2  = *.*3, 
and 

*,=  3aGH0K/2 [(60GHO  - a3/?3  + OtfAfj)  + io(6/lo  + fi0«J/Wo)l- 


) (44b) 


\ 


(45) 


} (45a) 


It  is  important  to  note  that  the  surprising  uniqueness  is  achieved  by  requiring  a regular  integral  at  the  North 
Pole  (Eq’s.  41  or  42)  and  by  truncating  the  series  of  the  solution  and  the  differential  equations  after  the  second 
power  of  sin  0.  Of  course,  without  the  second-order  truncation,  uniqueness  can  no  longer  exist.  Undetermined  co- 
efficients become  available  to  satisfy  prescribed  boundary  data,  say.  at  distant  continental  shorelines. 


] 


11 


In  the  present  global-tide  model,  the  arctic  tides  (EqY  43,  44,  and  45)  will  be  'onsidered  valid  up  to  1° 
colatitude  (the  first  land  occurs  at  colatitude  7°;  see  Bathymetric  Tables  in  Schwiderski  (1978a)).  For  colatitudes 
2°  and  3°,  a linear  interpolation  will  be  used  to  match  the  polar  solution  with  the  numerical  solution  computed 
south  of  3°  colatitude. 


It  is  interesting  to  observe  that  if  the  Coriolis  force  and  eddy  dissipation  are  neglected  (ft  = 0.  A0  = 0),  then 
the  second-order  arctic  tides  (Eq’s.  43,  44,  and  45)  become  exact  global  tides  with  the  constants 


K = aGH0Kl[(bl}GH0  - o2R2)  + ioB0RzIH0] 


K{  = 2K.  K2  = -2oRK,  *3  = 0. 


3 xr 


= -2k,  k2  = o.  rtj  = o. 


From  Equations  43,  44,  and  45,  one  concludes  that  all  second-order  arctic  tides,  f,  vanish  at  the  North  Pole  with 
the  same  order  as  their  corresponding  driving  equilibrium  tides  rj.  Hence,  only  long-period  tides  exist  at  the  North 
Pole. 


5.  THE  DISCRETE  OCEAN-TIDE  EQUATIONS  (DOTEs) 


A.  THE  1°  BY  1°  GRADED  GRID  SYSTEM 


With  the  exception  of  Antarctica  (south  of  colatitude  9 = 168°),  .ue  entire  (ocean  and  land)  area  of  the 
globe  is  covered  by  a 1°  by  1°  grid  system  that  is  “graded”  toward  the  poles.  Each  spherically  rectangular  mesh 
cell  Sm  n is  bounded  on  the  east  and  west  by  longitudes  \m  = m°  and  >m_  = (m  - n)°,  respectively,  and  to 

the  south  and  north  by  colatitudes  9n  = n°  and  9n_  j = (n  - 1)°,  respectively,  so  that 


where 

m = p,  2/i, . . . , (360  -K)) 
and 

n = 1,2 168. 

The  “grading”  parameter  n = nn  is  defined  by 

li  = 1 for  n = 30  to  150 
H = 2 for  n = 15  to  29  and  n = 151  to  168 
li  = 4 for  n = 8 to  14 
H = 8 for  n = 1 to  7 


(46) 


(47) 


As  was  mentioned  in  Section  2(f),  the  grading  of  the  network  toward  the  poles  is  necessary,  in  order  to  main- 
tain a more  uniform  mesh  area  for  higher  accuracy  and  stability  (Sect.  5.G)  of  the  discrete  tide  model.  In  fact,  the 
grading  equations  (Eq’s.  48)  have  been  chosen  in  such  a way  that 


/usinn°>^forn  = 4,8,  15.30.  150,  165;  (49) 

i.e.,  the  southern  mesh  size  remains  larger  than  half  the  equator  mesh  size.  This  important  condition  is  violated  for 
n = 1°,  2°,  3°  and  n - 166°,  167°,  168°.  However,  in  Section  4 it  was  pointed  out  that  the  numerically  discrete 
tide  model  begins  at  colatitude  9 = 4°.  For  colatitudes  9 = 3°  and  2°.  the  numerical  solution  will  be  matched  by 
linear  interpolation  to  the  second-order  arctic -tide  approximations  derived  in  Section  4.  which  are  assumed  to  be 
useful  up  to  colatitude  9 = 1°  As  will  be  shown  in  Section  5.G.  the  slight  violation  of  the  condition  in  Equation  49 
at  the  three  southern  colatitudes  is  not  so  severe  as  to  affect  the  stability  characteristics  of  the  model.  In  any  case, 
no  need  was  apparent  to  justify  an  additional  grading  step,  which  is  accompanied  by  an  unnecessary,  extra  compu- 
tational effort  (Sect.  5.D). 


In  the  global  network  defined  above,  land  and  ocean  mesh  cells  are  distinguished  by  zero  or  nonzero  depth 
data  which  will  be  assigned  to  each  cell  (Sect.  5.B).  The  “mathematical”  boundaries  of  the  oceans  follow  in  an 
obviously  zigzagging  fashion  the  mesh  lines  of  boundary  cells.  However,  empirical  tide  data  known  at  continental 
and  island  stations  will  be  utilized  to  specify  indirectly  more  physically  true  boundaries  (Sect's.  2(0  and  5.F). 
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PRECEDING  PADS  BLaMC-NOT  FILMED 


B.  HYDRODYNAMIC AL  OCEAN  BATHYMETRY 


The  ocean-depth  data  collected  by  S.  M.  Smith  et  al.  (1966)  were  rearranged  and  linearly  interpolated  to  fit 
the  new  1°  by  1°  graded  grid  system  described  in  Section  5.A.  The  original  data  bank  had  to  be  corrected  for  ob- 
vious errors  in  continental  and  oceanic  labeling,  ocean  and  land  signs,  shorelines,  and  some  exponents.  All  land 
elevations  were  set  to  zero  including  all  depth  data  less  than  50  m.  Furthermore,  the  following  meshwise-discon- 
nected  border  seas  were  excluded  from  consideration  by  assigning  zero-depth  values  to  their  corresponding  mesh 
cells:  the  Baltic,  Kattegat,  Irish,  Mediterranean,  Red,  Japan,  Sulu,  and  Ceram  Seas;  the  Hudson  and  Korean  Bays; 
and  the  Chihlian,  Persian,  and  Californian  Gulfs. 

After  some  negative  computational  experience,  the  depth  data,  which  were  originally  defined  by  applying  a 
''hydrostatical''  averaging  principle  (see  Sect.  2(0),  were  revised  using  the  following  “hydrodynamical”  principles: 

(a)  Boundary  cells  at  or  near  continental  shorelines  consisting  of  more  than  half  oceanic  areas  of  depths 
larger  than  5 m were  designated  ocean  cells,  and  the  average  oceanic  depth  values  were  assigned  as  the  “hydro- 
dynamically”  averaged  depths  of  the  entire  cells.  The  new  depth  value  is  preferable  to  the  “hydrostatically”  averaged 
depth,  which  preserves  the  actual  water  mass  but  ascribes  artificially  a shallow  shelf  character  to  the  cell. 

(b)  Island  cells  were  declared  terrestrial  cells  with  depths  zero  if  either  the  island  areas  were  larger  than  half 
the  mesh  areas  or  the  (elongated)  island  lengths  exceeded  the  mesh  diameters. 

(c)  Island  cells  that  remained  oceanic  cells  were  assigned  depth  values  less  than  the  hydrostatically  averaged 
values.  In  this  case  and  in  situations  of  submerged  seamounts  or  narrow  ocean  ridges  (e.g.,  Aleutian,  Marianas,  and 
Carribbean),  the  hydrodynamical  depths  depended  on  the  assessed  “barrier”  effects  of  the  current  obstacles:  the 
longer  and/or  higher  the  barrier,  the  lesser  the  depth. 

(d)  The  assigned  minimum  depth  was 


Hm  = min //(A,  0)  = 20  m,  (50) 

which  is  lowered  to  10m  by  the  averaging  equations  (Eq’s.  7 1 ).  The  maximum  depth  was  set  at 

Hm  = Max  //(X,  0)  = 7 000  m,  (51) 

which  eliminated  a few  isolated  deeper  values  that  unnecessarily  lessened  the  stability  of  the  tide  model  (Sef.  5.G). 
The  averaged  North  Pole  depth  was  found  to  be 

H0  = 3 600  m,  (52) 

which  is  needed  to  compute  second-order  arctic  tides  (Sect.  4)  for  colatitude  line  n = 1 . 

(e)  All  depth  data  H n are  considered  representative  for  the  center  of  the  cell  Sm  n ; i.e.,  for  m = n,  2#i,. . . , 
360  and  n = 1,2,...,  168, 

Hm.n  ="(W>  (53) 

where  \m  = (m  - p/2)°  and  0n  = (n  - 'ti)° . 
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As  was  already  pointed  out  in  Section  2(0.  the  hydrodynamically  justified  principles  discussed  above  in 
(a)-(e)  are,  naturally,  quite  subjective  and  by  no  means  free  of  any  error.  Nevertheless,  some  computational  experi- 
ments indicated  only  very  minor  effects  of  isolated  depth  data  changes.  More  than  3 000  depth  values  were  changed 
but  only  very  few  of  those  required  additional  readjustments  in  order  to  keep  some  limitation  on  the  first  and 
second  derivatives  of  H(\,  0);  i.e.,  on  the  relative  differences  given  by  Equations  71.  Furthermore,  the  hydrody- 
nantical  interpolation  of  empirical  tidal  data  (Sect.  5.F)  known  at  continental  and  island  stations  greatly  diminishes 
the  need  for  precise  boundary-depth  data.  The  revised  depth  data  bank  used  in  the  present  tidal  computations  will 
be  published  in  Schwiderski  (1978  a). 


C.  EMPIRICAL  TIDE  DATA 

The  present  tide  model  incorporates,  by  a unique  hydrodynamical  interpolation  procedure  (Sect.  5.F),  empiri- 
cal tidal  data  observed  and  harmonically  analyzed  at  numerous  continental  and  island  stations.  These  data  were 
taken  from  publications  by  the  National  Ocean  Survey  (1942),  the  International  Hydrographic  Bureau  (1966),  the 
British  Admiralty  (1977),  and  by  Pekeris  and  Accad  (1969),  Zahel  (1970,  1973),  Cartwright  (1971),  and  Luther 
and  Wunsch  (1974).  Unfortunately,  the  most  recent  publication  by  the  British  Admiralty  (1977)  lists  harmonic  con- 
stants only  for  the  four  major  tide  components  M2,  S2,  Kp  and  Oj  and  excludes  the  European  waters  completely. 

The  voluminous  data  banks  had  to  be  screened  in  order  to  eliminate  observations  that  are  meaningless  or 
unreliable  for  the  present  ocean-tide  investigations.  For  example,  tidal  constants  were  excluded  that  were  listed  for 
stations  deep  inside  estuaries  or  narrow  bays  (e.g.,  Hudson  River.  Bay  of  Fundy),  at  the  mouth  of  large  rivers  (e.g., 
Amazon),  between  sheltering  islands  (e.g.,  Alexander  Archipelago,  Solomon  Islands),  and  inside  sheltering  reefs 
(e.g..  Great  Barrier  Reef). 

About  2 500  stations  were  selected  for  further  examination  of  their  data  concerning  locally  restricted  dis- 
tortions. For  instance,  some  data  taken  over  short  distances  along  a coastline  displayed  rather  drastically  alternating 
times  of  high  water,  which  are  obviously  meaningless  for  oceanic  tidal  studies.  At  many  stations,  different  tables 
gave  different  tidal  constants.  Some  of  those  discrepancies  at  island  stations  are  shown  for  the  M2  tide  in  Table  2. 
Similarly,  for  some  mesh  cells,  several  different  station  data  were  available,  and  only  one  representative  average  had 
to  be  chosen.  This  situation  is  illustrated  in  Table  3 for  the  M2  tide  around  Bermuda.  Many  of  those  differences 
can  probably  be  explained  as  simple  errors  in  printing  or  computing.  For  instance,  the  phase  difference  of  about 
1 hr  at  Port  Galets  on  La  Reunion  Island  ( Table  2)  seems  to  be  due  to  some  error  in  observing  the  correct  reference 
time,  which  varies  from  listing  to  listing.  Most  differences,  such  as  those  shown  for  Bermuda  stations  in  Table  3, 
are  definitely  true  local  variations.  In  this  connection,  the  important  tidal  measurements  by  Gallagher  et  al.  (1971) 
at  Fanning  Atoll  in  the  Central  Pacific  may  be  mentioned.  Tides  outside  and  inside  the  small  atoll’s  lagoon  differed 
by  about  50%  (20  cm)  in  amplitude  and  by  a phase  lag  of  about  50°  (1  hr,  40  min.). 

In  general,  the  most  recent  listings  in  the  British  Admiralty  Tide  Tables  (1977)  were  chosen  over  older  tabula- 
tions as  most  reliable.  The  selection  of  the  data  was  further  aided  by  the  earlier  and  subsequent  tidal  computations. 
Altogether,  some  1 700  M,-tide  data  were  selected  and  assigned  to  the  centers  of  their  respective  mesh  cells.  Using 
linear  interpolation  and  tidal  computations,  the  total  number  of  prescribed  tide  data  used  in  the  M2-tide  construc- 
tion was  increased  to  more  than  2 000.  Essentially  all  continental  boundary  cells  carry  empirically  supported  tide 
data.  The  empirical  coverage  is  only  marginal  at  Arctic  and  Antarctic  shorelines.  Most  empirical  title  data  known  a' 
island  stations  are  also  included  in  the  tide  model.  All  empirical  tide  data  will  be  distinguished  from  computed 
values  by  underlining  in  the  printed  tide  tables,  which  will  be  published  (see  Tables  5-8)  as  supplemental  parts  of 
this  report  (Schwiderski,  1978  b). 
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Tabic  2.  Empirical  Mj-Tidc  Differences 


Station 

Latitude.  Longitude 


Tenerife. Canary  Island  (A) 
28°29'N  I6°I4'W 


Pott  Praia.  Cape  Verde  Island  (A) 
I4°S5'N  23°3I'W 


Port  Calets.  La  Reunion  Isla.  (I) 
20°55'S  55°  1 7'L 


Mawson.  Antarctica  (I) 

67°36'S  62°53'L 


Wilkes  Station,  Antarctica  (II 
66°15'S  110°3l'E 


Welles  Harbor.  Midway  Island ( P ) 
28°I2'N  1 1 7°22’W 


Lniwetok  Atoll,  Marshall  Is.  (P> 
I l°2l'N  I62°2 1 'E 


lies  Wallis.  Fiji  Island  (P) 

I3°22’S  I76°H'W 


Suva  Harbor,  Viti  Levu.  Fiji  l.(P) 
I8°0^'S  I78°26'F 


B.A.T.  (77 )a  N O S,  (42)d  Others  Initialed 

b(  tit ) 6C(°)  £(  in ) 61°)  £(nt)  6(°) 


0.53 

178 


0.5b  0.50 

195  212 


‘‘B.A.T.  (77)  = British  Admiralty  Tables  ( 1977). 
= tidal  amplitude. 

= tidal  phase  relative  to  Greenwich. 
dN.O.S.  (42)  = National  Ocean  Survey  (1942). 
CZ  = Zahel  (1970). 

'P  = Pekeris  and  Accad  (1969). 
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Table  3.  Bermuda  M^-Tidc  Observations 

Station 


latitude.  Longitude 

f (cm) 

6b(°) 

Reference 

St.  George's  Island 

32.55N.  64.70W 

36 

359 

Bi  iltsh  Admiralty  (1977) 

St.  David's  Island 
32.37N.b4.65W 

— 34 

355 

Bntisli  Admiralty  ( 1977) 

Great  Sound 

32.32N.64.83W 

38 

6 

Br Uish  Admiralty  ( 1977) 

St.  George's  Island 
32.37N.64.70W 

35 

0 

U.S  Nat.  Ocean  Survey  ( 1942) 

St.  George’s  Island 

32.40N.  64.70W 

37 

0 

Pekeris  & Accad  ( 1969) 

St.  George's  Island 

32.37N,  64.70W 

36 

359 

Zahcl  ( 1970) 

St.  George's  Island 

32.40N,  64.70W 

36 

358 

/ettler  el  al  1 1975) 

Deep  Sea  (GOBI  IV) 

32.29N,  64.50W 

38 

1 1 

J T.  Kuo  Letter  ( 1977) 

= tidal  amplitude. 

b6  = tidal  phase  relative  to  Greenwich. 


Table  4a.  Deep-Sea  M->-Tidc  Data  lor  the  Gulf  of  Mexico  and  Caribbean  Sea 


Station 

Obseived  1 

1 Model 1 

1 Fttot 

Latitude.  Longitude 

(‘(cm) 

fiV) 

Stem) 

fil”) 

A$(cm) 

mC) 

W Florida  Shelf  St 

26.71N,  84.25W 

7 

97 

7 

92 

0 

5 

Deep  Gulf  St 

24.77N,  89.65W 

1.3 

226 

1.6 

225 

+0  3 

1 

Mistcriosa  Bank 

I8.88N,  83.81  W 

8 

84 

9 

89 

♦ 1 

♦5 

Rosalind  Bank 

16.61  N.80.34W 

m 

8 

102 

♦ 1 

5 

l-last  Carib.  St.  (6-month) 

16.54N,  64.88W 

urn 

n 

♦ 1 

-5 

Hast  Carib.  St.  (1 -month) 

I6.52N.  64.91  W 

0.6 

153 

1.5 

148 

0.9 

_5 

= tidal  amplitude. 

b<5  = tidal  phase  relative  to  Greenwich. 


J 
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Table  4b.  Deep-Sea  Mj-Tide  Data  for  the  Pacific  and  Atlantic  Oceans 


Station 

Latitude,  Longitude 

Observed 

Model 

Error 

>(cm)  «D(°) 

{(cm)  «(°) 

A£(cm)  A6(°) 

Pacific  St.  1 (Middleton) 

58.76N,  145.7  IW 

no 

284 

included 

Pacific  St.  3 (Tofino) 

48.97N,  127.29W 

99 

239 

included 

Pacific  St.  (San  Francisco) 

38.I6N,  124.91W 

54 

227 

included 

Pacific  St.  (Josie  II) 

34.00N,  144.99W 

27 

267 

27 

273 

0 

+6 

Pacific  St.  (Flicki) 

32.24N,  120.86W 

43 

149 

included 

Pacific  St.  (Josie  1) 

31.03N,  I19.80W 

43 

142 

included 

Pacific  St.  (Kathy) 

27.75N.  124.37W 

29 

128 

27 

130 

2 

+2 

Pacific  St.  (Filloux) 

24.78N.  I29.02W 

19 

107 

18 

105 

-1 

_ *> 

Atlantic  St  1 (N  Y Bight) 

39.32N,  64.36W 

44 

350 

included 

Atlantic  St.  (N.C.  St.  1 ) 

32.69N.  75.62W 

48 

356 

46 

358 

■> 

+2 

Atlantic  St.  (Savannah  B) 

31.95N.80.68W 

88 

15 

included 

Atlantic  St.  (Scope) 

30.43N,  76.42W 

45 

358 

46 

3 

+2 

+5 

Atlantic  St.  (AOML  1 ) 

28.14N,  69.75W 

34 

1 

35 

6 

+ 1 

+5 

Atlantic  St.  (AOML  3) 

28.24N,  67.54W 

34 

359 

34 

4 

0 

+5 

Atlantic  St.  (MFRT) 

27.99N,  69.67W 

34 

360 

34 

6 

0 

+6 

Atlantic  St.  (REIKO) 

27.97N,  69.67W 

35 

1 

34 

6 

1 

♦5 

Atlantic  St.  (EDIE-May) 

26.46N.  69.33W 

32 

3 

32 

7 

0 

+4 

Atlantic  St.  (EDIE-March) 

26.45N,  69.32W 

31 

1 

32 

7 

+ 1 

♦6 

= tidal  amplitude. 

b6  = tidal  phase  relative  to  Greenwich. 
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Naturally,  it  must  be  remembered  that  the  selection  of  representative  empirical  tidal  data  (compare  depth 
data.  Sect.  5.B)  is  not  at  all  free  of  subjective  judgment  and  may  be  somewhat  erratic.  Obviously,  only  future  ad- 
ditional tidal  measurements  can  improve  this  model  in  this  respect.  Nevertheless,  according  to  the  instruction  notes 
accompanying  the  British  Admiralty  Tide  Tables  (1977),  it  can  probably  be  assumed  that  almost  all  important  tide 
data  selected  carry  an  accuracy  that  is  at  least  as  high  as  the  desired  10  cm  specified  in  Section  1.  In  any  case, 
computational  experiments  showed  that  isolated  reasonable  variations  of  the  boundary-tide  data  do  not  affect 
significantly  the  adjacent  oceanic  tides.  It  was  also  found  insignificant  to  the  overall  quality  of  the  tide  model 
whether  the  empirical  data  were  assigned  to  the  centers  or  to  the  shore  boundaries  of  the  respective  cells  (see 
Part  II). 

Attempts  were  made  to  incorporate  also  recent  deep-sea  tidal  measurements  (Sect.  2)  into  the  present  model. 
Since  the  hydrodynamical  interpolation  of  empirical  data  is  essentially  based  on  bottom  and  boundary  irregularities 
(see  Sect.  5.F(a)-(d).  (a)),  no  physically  valid  justification  was  found  to  include  distant  offshore  deep-sea  measure- 
ments into  the  model.  However,  some  deep-sea  measurements  near  rough  shore  and  bottom  areas  were  included. 
Fortunately,  without  exception,  all  excluded  offshore  deep-sea  measurements  known  to  the  author  agree  very  well 
with  the  computed  M2-tide  data  (see  Table  4). 

D.  DERIVATION  OF  DISCRETE  OCEAN-TIDE  EQUATIONS 

Following  essentially  the  hydrodynamical  numerical  method  of  Hansen  (1966)  and  Zahel  (1970),  the  COTEs 
(Eq’s.  29  and  30)  may  now  be  converted  to  an  explicit  finite-difference  analog  called  “discrete  ocean-tide  equations" 
(DOTEs).  In  a first  step,  all  spacial  derivatives  will  be  replaced  by  divided  central  finite  differences  using  the  Richard- 
son (1922)  staggered  scheme  (Sect.  2(g))  illustrated  in  Figure  2c.  In  agreement  with  the  graded  grid  system  defined 
in  Section  5. A,  it  is  convenient  to  introduce  the  following  notations: 


A 0 = jt/360  = 1/2°  mesh  size. 

(54) 

A X = /r  A 0 , L = 2RA9, 

(55) 

A t = time  step  to  be  specified; 

and  for  a fixed,  oceanic  mesh  cell  Sm  n (Eq’s.  46,  47,  48)  at  u-points  (Figure  2c) 

\n(“)  = 2(m  - #i)AX , 

6„(u)  = (2 n - 1)A0  , 

’»£.“)  = »?(*,*(“)’ M“).  0. 

r„(“)  = l//a  sin  0„(m); 


(56) 


at  v-points  (Figure  2c) 


\„(v)  = (2m  - *i)AX, 

6 (v)  = 2wA0. 

no  = »'(\„(vM„(v).f). 

m.n  " 

v (t;  v)  = n (Xm(v),  0„(v),  /). 
rn(v),=  l/#i  sin  0„(v); 
and  at  f-points  (Figure  2c) 

= Vv>- 

0n<f>  - 

m * t(xM(v).  «„(«).  o 

m.n 

For  the  first  COTE  (Eq.  29a)  at  u-points  one  has  the  following  differential-difference  replacements. 


(57) 


(58) 


4AX  U (/)  - U (0  - V (t) 

X.m.n  m+n.n  m - p.n 


4AXJ  U (/)  -*  V (/)  - W(t)  + U ( t ) 

KX.m.n  m + n,ii  m.n  m-n.n. 


4 AOU  (/)  - U (/)  - U (/) 

O.m.n  m.n*  1 m.n- 1 


4Ad2U  (I)  -*  U (/)  -2 (/(/)  +(/(0 

06. m.n  m.n  + 1 m.n  m.n-  1 


4^  (r;u)  - [F  (/)  + K(r)  1 + [^  (/)  + K (0 

m.n  m.n  m.n- 1 m-ji.n  m-n.  n- 1 


4AX  F (t.v)  - [('(Z)  + K(/)  | - |F(/)  + ^ (r) 

X.m.n  m.n  m.n-  1 m-/i,n  m-n.  n- 1 


2AXf(f.u)  - f (/)  - f (/) 
X.m.n  m.n  m-n.n 


(59) 
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For  the  second  COTE  (Eq.  29b)  at  v-points  one  has  the  analogous  differential-difference  substitutions 
4AAJF(r)  -*■  V (/)  - 2V  (0  + V (r)  , 

W.m.n  m+n.n  m.n  m-n,n 


4AdV  (0  -*V(t)  - V(t) 

8, m.n  m.n*  1 m.n- I 


4A62V  ( t ) -*  V (t)  - 2 V(0  + V (t)  , 

00.  m.n  m.n  + l m.n  rn.n-1 


4V  (t;  v)  - [0  (/)  + 0(0  1 + [0  (l)  + 0(0  J , 

m.n  m+n.n  + 1 m.n  m+n.n  m.n*  1 


4A\0  (/."v)  - \U  (0  ~ 0(0  ) + 1 0(0  - 0 ] , 

K.m.n  m+n.n+l  m.n  m+n.n  m.n*  1 


4A60  (t.y)  -*  [0(t ) - 0(0  \ - [0(0  - 0 [. 

8. m.n  m+n.n  + l m.n  m+n.n  m.n+l 


2A0f  (r.v)  -*  HO  - t (0  • 

O.m.n  m%n+ 1 m.n  ) 

The  third  COTE  (Eq.  30)  at  f -points  becomes  directly 

R sin  0 ,(n)  f (/)  =-±-[0(0  -0(0  ] + [sin  8MV(t)  - sin  9 (v)  V (/)  ].  (61) 

" t. m.n  -A*  m.n  m+u,n  m.n  n-  1 m.n- 1 

It  may  be  noted  that  in  Equations  59  and  60  V(t;  u)  at  u-points  and  0(t;  v)  at  v-points, are  defined  as  linear 
averages  of  adjacent  V(t)  at  v-points  and  0(0  at  w-points.  respectively.  After  defining  obvious  analogous  substitu- 
tions for  the  depth  function,  H(\.  6)  (given  at  f-points.  Sect.  5.B),  and  its  first  and  second  partial  derivatives  at  m- 
and  v-points,  one  arrives  at  the  following  spacially  differenced  equations: 

GH  ( u ) 

AtO  (0  = aAt  ” n(t.u)  + A{i)  [?  (f)  - f (0 1 ~ 

t.m.n  a s,n  Vn\U)  K.m.n  m.n  m-n.n  m.n  m.n  m.n 


A5 

0 (0 

+ A 6 0(0 

+ A 7 

o(o  + ■ 

A*0  (r) 

m.n 

m+n.n 

m.n  m-n , 

7i  m.n 

m.n*  1 

m.n 

A9 

l(V(0 

+ v (r)  , ) - 

(V  (0 

+ v(t) 

)] 

m.n 

m.n 

m.n-  1 

m - n.n 

m- n.n- 

i 

A 10 

uno 

+ v (0  ) + 

(V(0 

+ v(0 

)). 

m.n 

m.n 

m.n-  1 

m - p.n 

m - n.n- 

- 1 

I 


GH  (v) 

AfK(0  =<* At-^-  rj(/.v)  + B3  |f (/)  - f(r)  ] - B*  V(r) 


t.m.n 


8. m.n  m.n  m,n  + 1 m.n 


+ B5  [V  ( t ) + F (r)  1 + B6  K(r)  + B7  F(r) 

m.n  m*n.n  m- n.n  m.n  m,n+l  m.n  m.n- 1 

+ B*JU  (/)  - U(t)  | + fl’  |t/(/)  -{/(/)  | 

* m+ji.w+l  m.n  w+/j,/i  m.n  + l 


+ (0  + U(t)  + (/(f)  + (/(/)  ] , 

• m+u.  h+1  m.n  m +ji,n  m.n  + 1 


<6  2b) 


and 


Af?  (f)  = C‘  [(/(f) 

r.m.rt  m,/i 


€/(r>  1 + C^(0 

m+n,n  m.n 


cy«) 


m.n  I 


(63) 


where  all  A.  B.C,  and  H's  are  defined  by  Equations  69,70,  and  7 1 . 

In  the  second  part  of  the  differencing  process,  a regular,  finite-difference  scheme  (Sect.  2(g),  Figure  3a)  in 
time  is  defined  by  integrating  Equations  62  and  63  over  a single  time  step.  At,  using  an  average  integration  rule  of 
the  form  ((/  -*  V-*t) 


where 


f'/+1 

I U(t)dt  --  At 


k(/'+1  +(1  -K)U' 


(64) 


u'  = mtj). 

tj  = (j  - l)At,i  =1.2 (65) 

and  k - some  differencing  parameters,  which  usually  satisfy  the  restriction 

0<x<l.  (66) 


As  will  be  seen,  the  resulting  discrete  ocean-tide  equations  (DOTEs)  are  very  sensitive  to  the  choice  of  the  differ- 
encing parameters,  x,  and  the  time  step.  At.  In  fact,  depending  on  the  chosen  values  of  k and  At,  the  DOTEs  may 
be  stable  or  unstable  (Sect.  5.G)  for  any  specified  values  of  eddy  viscosity,  A,  and  bottom-friction  coefficient,  B. 
Moreover,  different  values  of  k for  Equations  62a,  62b,  and  63  and  for  the  various  point  values  of  U,  V,  and  J can 
be  chosen  so  that  the  resulting  DOTEs  may  become  explicit  or  implicit.  Considering  the  complexity  of  the  ocean 
basin  and  the  large  number  of  oceanic  mesh  cells,  an  explicit  form  for  the  DOTEs  is  here  preferred. 
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Alter  carrying  out  the  integration  of  Equations  62  and  63  with  the  three  (obvious)  differencing  parameters 
k = 0.  k . and  if.  one  arrives  at  the  following  explicit  DOTEs: 


l1  + Um.n  = K,.nsin2T^  ~ *>  + Am.n  cos^(2/'  - 1) 

+ k,u„  - + ii  - (i  - 


+ As  U>  + A6  U 1 + A1  U'  + A6  U' 

m.n  m+n.n  m.n  m-v.n  m.n  m,n  + \ m.n  m.n-  1 


+ + C»-l>  " + VL>.n- 1>1 


+ + VL.n-{>  + lVL*.n  + 

['  + cos£f(2/'  - ')  + ®m.n  s*n~y—  (2/  - 1) 

+ - o + ii  - (i  - *<„) 

+ + VL^  + + , + Bm.nVi.n-X 


(67a) 


and 


+ fl8  \uf  - U'  1 + fi9  \U'  - U'  .1 

m,n  1 itv +j«(ff  + l m.n1  m.n*  m+p.n  m.n  + \ 1 


+ fl10  ((/'  + u1  + u'  * U1  . 

m.n1  m*n.n+ 1 m.n  m*n.n  m.n*  1* 


C = t.n  **KK'n  - O + r»C!  - 


(67b) 


+ (I  -«)IC,K.W  - <£♦„.„>  * 


(68) 


The  coefficients  of  the  DOTEs  (Eq’s.  67  and  68)  are 


Al.„  = ^m.«COs2*'(m  - /i)A0’ 

Al,n  = Am.nSin2v(n  ~ ti)A6’ 

A]nn  = P&tGH  (u)  r (u)/L, 

m.n  mn 

Am.n  - 2*f  ' "(“)!  + • 

^ ffi.w  wi./i 

4*  11  + 

m,n  L n n mn  m.n 

4m»  *aT'W“)r>)"<“>  l1  -"(«>]• 

’ **  m,#t  m.n 

47  n = flT  H (M)  1<Mm>  + ^ <2  + r (M))  COS  (2ft  - 1)A0|  , 

• l-  m.n  L 

4?  ~ aT  H(U)  ^ (2  + r (m))  COS  (2/1  - 1>A0]  , 

• l.  nt,n  *■ 

Al,  „ = a^r<u)//(M)[t//  <«)//<«)  - ^(3  + 2r(u))  cos  (2/*  - 1)A0], 
L m.n  m.n  * 

and 

Amn  = -«T  r<M>w<“)  "<">  (*„<«)//(«)  -^C0S(2 //  - U A0| 

• n m.n  m.n  m.n 

+ ^ J2  cos  ( 2/t  - DAO  , 


where 


(69a) 


(69a) 


Bl,n  = Bn,,,C0%V(2m  ~ »)A0- 

Bl.„  = -m)A6. 

fl3„  „ = P&tGH (v)  IL. 

m,n  m.n 

B*n  = 2a  ^ !i/  (v)[w  (v)//(v)  - tf(v)|  ♦ 6Af/r(,(v)//(v)  . 

• m.n  m.n  m.n 

=a  X*.(v>r2<v)WjV>  ’ 

, = «7ff<v)l*,(v)(l  - W(v))  + ^(3  + 2r  (v))cos2/»Afl|, 

* l*  m.n  m.n  * 

Bl  „ = [*  (v)(l  + W(V)  ) - ^ (3  + 2r  (v))cos2«A0]  , 

in.  #f  1 m.w  ^ 

Sm  « = aT  UnAflr  (v)cos2nA0  - tf(v)  ] . 

Bmn  = afr»WMv>w(v)  [2/iAflr  (v)cos  2nA0  + /7(v))  , 

’ m n ' m . n 


Bn?„  = rn(v>*»M//(v)WW  ItiAOr  (v)cos  2/zA0  ♦ tf  (v)  1 

•-  m.n  m.n  m.n 


- — £2  cos  2«A 0 


where 


m ,n 


.ccG  . . . a At 
4 Zb  A://(v)sm-r- 
OK  m.n  l 


2sin4nA0,  v - 2 

cos  4nA0,  v * 1 

-3  sin  4nA0,  v = 0 


and 


+„(*)  =3(1  + l/r„(v)),  w„(v)  = 1 t r>  (v); 


(6%) 


(69b) 
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<vi  =4r..(,,)- 


- /i-y-r„(«)  Sin  2 nAO. 


Cl  = |U-j-T)((tt)  sin  2(«  l)A0. 


In  these  coefficients,  the  following'depth  functions  were  used  at  w-points: 


"}">  =Tl/V»  + I- 


//(,<)  = -#(«>  1. 

m.n  m*n.n  m-fA.n 

m.n 


H{,,)  = -T77777T  l7/  <“>  ~ //<">  1 • 

m.n  in.it—  1 »H.ii+l 

in  ,ii 


//  (M)  = H (M)  iw  (M)  2ri 2(w)  [//  (w ) + (nAd)2 

m.n  m.n  ( * ' m.n 


+ /iA0^  (u)  cos  (2;i  - 1)A0  [(3  + 2r  (m))//(w)  - pA0T  (w)  cos  (2/i  - l)A0)i 

iii.ii  ; 


-y[r2(«)(//  (n)  +//<«)  ) + (//(«)  + //  (m)  )]; 

- in+fi.ii  in -jj.n  m.n*  1 ni.w-l 


and  at  v-points 


"<v„>  =Tl//».»  + 'Vn+iJ  • 


/7(v)  =17/1^  |//(v)  - H (v)  1 • 

m.n  *+ri  \v p m*n.n  m-pi.n 

m . n 


//(v)=  4/77^|/y  (v)  , W<V>1- 

iii.ii  \v  p m.n- I m.n*  1 

m . n 


W(V)=  //(V)  L (v)  2|//*(v)  + (/iA0)2 r2(v)l 

m.n  m.n  t m.n  " 

+ f/AOiJ/  (v)/‘/(v)  [3  + P (v)]  cos  2iiAd\ 
' m.n  " ' 


i |r2(v)(W(v)  + H (V)  ) + [H  (v)  + H (v)  )1. 

“ iii*ii.ii  in-iA.ii  iii.ii*  1 m.n- I 
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It  may  be  mentioned  that  for  x =0  and  5T  = I , the  finite-difference  scheme  coincides  with  the  technique  used, 
e.g.,  by  Hansen  (1966),  Zahel  (1970,  1973,  1975),  and  Estes  (1975,  1977).  Extensive  exploratory  computations 
were  carried  out  by  the  author  with  numerous  k and  if  values  within  the  ranges  0 < x <1.5  and  0.5  < x < 1 . The 
computations  produced  no  drastic  differences,  provided  the  eddy  and  bottom-friction  coefficients  a and  b and  the 
time  step  Af  were  suitably  chosen  within  their  respective  stability  constraints  (Sect.  5.G).  Finally,  it  was  decided 
to  use  the  values 

k = 7T  = 1 (72) 

because  they  yielded  a preferable  stability  and  seemingly  best  results.  Moreover,  this  deviation  from  the  Hansen- 
Zahel  method  became  most  significant  for  the  hydrodynamical  interpolation  of  empirical  tidal  data  described  in 
Section  5.F.  Indeed,  this  novel  technique  uses  the  special  property  that  for  x = 1 the  bottom-friction  coefficient, 
b.  which  enters  only  in  A 4 and  B 4 (Eq’s.  69),  becomes  essentially  a scaling  multiplier  of  and  F^1,,  in  Equa- 
tions 67a  and  67b.  Thus,  together  with  if  = 1 in  Equation  68.  the  bottom-friction  coefficient  can  easily  be  adjusted 
locally  to  match  more  closely  prescribed  tidal  data. 

With  x and  x"  specified  by  Equation  72,  the  DOTEs  (Eq's.  67  and  68)  still  contain  the  parameters  a,  b,  and 
A 1.  which  remain  at  one's  disposal  within  their  respective  stability  ranges  (Sect.  5.G).  They  will  be  utilized  to  achieve 
best  results  by  trial-and-error  computations.  The  DOTEs  (Eq's.  67a,  b,  and  68)  can  be  applied  to  all  oceanic  mesh 

cells  Sw  n with  m = p,  2^t 360  and  n = 4,  5 168  sweeping  across  the  globe  from  n = 4 to  n = 168.  This 

procedure  can  be  executed,  provided  suitable  initial  and  lateral  boundary  data  (Sect.  5.E)  are  prescribed.  At  co- 
latitude line  n = 4,  the  numerical  solution  is  matched  to  the  second-order  arctic  solution  (Sect.  4)  by  the  linear 
interpolation  (m  = p.  2/j 360) 

U'*\  = 1 PI//*1  + U'*x  | 
ni  .3  3 * *"  Hw . 4 rn , P * 

///♦  i = * i rji*  i + ->r//+ 1 1 

Sir. 2 3 'Sir. 4 -Sn.l1' 

where  Ul*\.  , . and  f£‘  are  computed  by  Equations  43,  44.  or  45.  For  colatitude  lines  n = 7,  14.  29,  1 50, 

spacially  corresponding  data  [V.  V,  f)  on  n =8,  15,  30,  151  (see  Eq's.  46,  47,  and  48)  are  defined  by  linear  inter- 
polation. Vice  versa  for  n = 8.  15,  30.  151 . spacially  corresponding  values  ((/,  V,  f)  on  n = 7,  14,  29,  1 50  are  also 
defined  by  linear  interpolation. 


E.  LATERAL-BOUNDARY,  INITIAL.  AND  FINAL  DATA 

In  order  to  complete  the  ocean-tide  model,  the  DOTEs  (Eq's.  67  and  68)  must  be  supplemented  by  suitable 
lateral-boundary  and  initial  values.  As  was  explained  in  Section  2(d),  in  turbulent-flow  situations,  the  mathematical 
boundary  conditions  usually  preferred  are  (a)  no-flow  across  the  ocean  shorelines  and  (b)  free-slip  along  the  ocean 
shorelines.  It  is  clearly  at  this  point  that  the  great  attractiveness  of  Richardson's  ( 1922)  staggered,  finite-difference 
scheme  (Sect's.  2(g)  and  5.D)  manifests  itself  in  the  practical  simplicity  with  which  the  no-flow  and  free-slip  (or  no- 
slip) boundary  conditions  can  be  worked  into  the  model.  In  fact,  if  5 is  an  oceanic  boundary  cell,  then,  by 

- 


(t/-F-f) 


(73) 
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definition  (Sect.  5. A and  Figure  2c),  the  zig-zagging  mathematical  boundary  lines  follow  only  mesh  lines  and  pass 
only  through  velocity  ( u and/or  v)  points.  The  no-flow  condition  (a)  is  implemented  by  declaring  at  those  points 


U'+V  = 0 and/or  F>+1„  = 0, 
m,  n nt,  n 


(74) 


respectively.  If,  say,  the  v-point  of  Sm  n is  oceanic  and  the  v-pointsofSm_^  n and/or  Sm+^  n are  terrestrial,  then 
the  free-slip  condition  (b)  is  satisfied  by  reflectively  setting* 


yi* 1 = + yi*  l 

m - n,n  m.n 


and/or 

m+n,  n 


= +F'+1  , 
m.n' 


(75) 


respectively. 


The  mathematical  boundary  conditions  (a)  and  (b)  were  applied  by  Hansen  (1966),  Zahel  (1970, 1973, 1975), 
Estes  (1975,  1977),  and  others.  In  the  present  model,  these  conditions  are  strictly  enforced  only  at  boundary  points 
where  no  known  empirical  tidal  data  are  available.  However,  at  most  oceanic  boundary  cells,  empirical  tidal  con- 
stants are  defined  in  Section  5.C.  These  data  are  worked  into  the  present  tide  model  by  a unique  hydrodynamical 
interpolation  procedure  (Sect.  5.F)  that  requires  a controlled  violation  of  the  no-flow  condition  (a). 


The  construction  of  the  tide  was  started  at  / = 1 (^  = 0)  with  an  ocean  completely  at  rest ; i.e.,  with  the 
initial  values  (m  = y,  2y,  * . . , 360;n  = 1,2,...,  168) 


{/I  =y  1 = f l = o. 

m,n  m,n  *m,n 


(76) 


The  computations  were  carried  over  a prescribed  number  of  time  steps,  / = J (mostly  a quarter  period),  and  then 
printed  for  inspection  of  the  results.  With  or  without  program  and/or  parameter  changes,  the  computations  were 
restarted  using  the  latest  or  any  earlier  taped  output  instead  of  the  initial  values  (Eq.  76).  Occasionally,  it  was 
beneficial  to  speed  up  an  unwanted  slow  decay  of  transient  eigenmodes  by  “negatively”  averaging  the  output  data 
of  half  a period  time  difference;  i.e.,  by  setting 


K,.n  = *[v,m,n-V’m:Tn (77) 

where  oAtJ  = it.  This  simple  procedure  diminishes  all  undesirable  eigenmodes  of  lower  frequency  than  the  forced 
frequency,  a.  and,  similarly,  most  higher  frequency  modes.  The  negatively  averaged  data  (Eq.  77)  represent  ob- 
viously improved  initial  data. 


Two  output  tidal  elevations  27^time  steps  apart  (mostly  a quarter-period  apart), 


ti.n  = COS  laAtJ  ~ VJ 

= «m,„  cos  [cAr(./~  27)  6m  J, 


(78) 


*lf  the  no-itip  condition  is  imposed,  then  the  + signs  in  Equations  75  must  be  changed  to  - signs. 
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were  used  to  compute  the  tidal  “amplitudes” 


««.«  = l*i  + 


2iV4 


(79) 


and  “phases” 


5m  n = arc  tan  y/x  ( = 0 for  x =y  = 0)  , 


(80) 


where 


and 


*1  = <nJ  + ^,„)/2cOSOA/Z 

= <«27-  fl„)/2sinaA/7 


y = f j sin  aAr(7  - /)  - £ 2 cos  gA/(/  - -f), 


x = £ j cos  oAt  (J  - J)  + ? j s'n  (/  ~ ■ 


(81) 


(82) 


In  the  finished  product,  i.e..  when  all  (unforced)  transient  eigenmodes  have  satisfactorily  decayed,  the  amplitudes 
(Eq.  79)  and  phases  (Eq.  80)  become  essentially  independent  of  time  and  undergo  no  further  variation  of  signifi- 
cance after  continued  computations. 

In  the  final  M2-tide  model  presented  in  Part  11,  the  convergence  toward  the  steady  state  of  the  amplitudes 
and  phases  was  found  to  be  generally  oscillating.  So  the  integration  process  could  safely  be  terminated  when  the 
amplitudes  and  phases  over  most  ocean  areas  varied  by  less  than  1 cm  and  1°,  respectively.  This  excellent  conver- 
gence feature  held  true  even  at  most  coastal  points.  The  convergence  was  slightly  less  complete  only  in  some  coastal 
cells  where  the  tidal  elevations  are  extremely  large  and  vary  rapidly  from  degree  to  degree.  But  even  in  those  anoma- 
lous places  (e.g.,  on  the  Patagonian  Shelf;  see  the  Tide  Tables  in  Part  11),  the  convergence  error  was  well  below  the 
desired  1 0 cm  specified  in  Section  1 . 

In  order  to  follow  the  convergence  of  the  computer  program  more  closely,  the  squared  tidal-amplitude  sum 


360 
■ 2 
m= 


H ^ m.n 


(83) 


was  computed  and  printed  for  each  fixed  colatitude  line  n - 4,  5,  ...  , 168.  To  compute  the  amplitudes  £m  n in 
Equation  83,  Equations  79  and  81  were  used  for  27=  1 and  J = /'  + I ; i.e.,  for  two  consecutive  time  steps  in  connec- 
tion with  Equation  68.  In  this  measure  (Eq.  83),  the  convergence  was  carried  to  almost  three  significant  figures 
for  all  n. 
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F.  hydrodynamical  interpolation  of  empirical  tide  data 


» 


The  discrete  ocean-tide  model  developed  so  far  is  completely  specified  with  the  exception  of  the  eddy  and 
bottom-friction  coefficients,  a and  b,  and  the  time  step,  At  (Eq’s.  69  and  70).  When  those  parameters  are  chosen 
uniformly  for  the  entire  ocean  basin  within  certain  stability  limits  (Sect.  5.G),  then  a unique  solution  can  be  com- 
puted, and  any  agreement  or  disagreement  with  observed  tidal  data  must  be  accepted  unless  a compromising  modifi- 
cation of  the  well-posed  problem  is  feasible. 


Preliminary  computations  carried  out  with  this  purely  theoretical  tide  model  produced  results  (Sect.  2(0) 
that  fulfilled  the  accuracy  requirements  posed  in  Section  1 over  most  ocean  areas.  Nevertheless,  a brief  glance  at 
the  many  significant  variations  of  ocean  tides  observed  at  continental  and  island  stations  (see,  e.g.,  the  underlined 
empirical  M^-tide  data  in  Tables  5-8,  and  Part  II)  leaves  oi.e  convinced  that  such  local  tidal  features  could  not  be 
accommodated  by  a mathematical  model  that  ignores  local  boundary  details. 


In  order  to  enhance  the  quality  of  the  purely  hydrodynamical  tide  model  described  above.it  appeared  desir- 
able to  seek  a physically  acceptable  way  to  incorporate  the  empirical  tidal  data  themselves  into  the  program.  In  a 
simple  exploratory  experiment,  the  arranged  empirical  tidal  amplitudes,  Xm  and  phases,  n (Sect.  5.C),  were 
introduced  into  the  computer  program  and  used  (when  available)  in  place  of  Equation  68  by  setting 


?/+1 

*m,n 


= ^m.ncos(oAti-bmn). 


(84) 


The  results  represented  a rather  spectacular  detailed  improvement  of  the  tidal  model,  very  close  to  the  final  tide 
data  presented  in  Tables  5-8  and  Part  II. 


Unfortunately,  the  substitution  of  Equation  84  for  the  sensitive  equation  of  conservation  of  mass  (Eq.  68) 
constitutes  a “continuity  gap” 


/+1 

m,n 


?'+l 

* m.n 


- f/+1 

*m,n 


(85) 


which  remained  bounded  but  grew  unrealistically  large  in  some  areas  after  continued  integration  in  time.  Also,  the 
velocity  field  near  such  observed  points  became  unacceptably  distorted. 


Of  course,  the  continuity  gap  (Eq.  85)  can  be  attributed  to  the  following  major  causes  which  are  physically 
plausible  and  have  been  briefly  mentioned  before: 

(a)  The  bottom-friction  coefficient,  b (in  A4  and  B 4 of  Eq’s.  69),  which  is  most  effective  in  boundary  cells, 
depends  on  local  shore  features  such  as  true  cell  size  and  bottom  slope  and  roughness  (Sect.  3.B(k)-(n)). 


(b)  The  boundary  cells  are  idealized  by  definition  of  strictly  mathematical  boundaries  (see  Sect's.  2(0  and 
5. A.  and  Figure  1 ). 


(c)  The  depth  data  of  boundary  cells  are  subjectively  defined  and.  hence,  faulty  (Sect's.  2(0  and  5.B). 

(d)  The  empirical  tidal  constants  in  Equation  84  are  also  faulty  to  some  degree  because  of  inaccurate  measure- 
ments. harmonic  analyses,  and  subjective  selections  and  assignments  to  the  centers  of  the  boundary  cells  (Sect.  5.C). 

(e)  The  discrete  ocean-tide  model  is  certainly  not  an  exact  description  of  the  true  oceanic  tide;  e.g.,  at  boun- 
daries, nonlinear  inertial  terms  assume  significance. 
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Obviously,  the  last  two  (hopefully  minor)  faults  can  be  reduced  only  through  continued  future  observations 
and  modeling.  However,  the  first  two  faults,  (a)  and  (b),  can  be  weakened  by  "hydrodynamically  interpolating" 
the  empirical  tidal  elevations  (Eq.  84)  into  the  tidal  model  and  narrowing  the  continuity  gap  (Eq.  85)  to  an  accept- 
able level  as  follows: 

(a)  Adjusting  the  velocity  field  by  a locally  controlled  implicit  variation  of  the  bottom-lriction  coefficient, 
b.  in  Equations  69. 


(E)  Lifting  the  strict  condition  of  no-flow  across  the  mathematical  ocean  boundary  and  allowing  for  a moni- 
tored in-  or  out-flow  by  implicitly  defining  a more  physical  ocean  boundary  (Figure  1 ). 

Since  the  finite-differencing  parameters  k and  kin  Equations  67  and  68  have  been  chosen  by  Equation  72, 
the  bottom-friction  coefficient,  b,  in  A 4 and  B 4 of  Equations  69  can  be  considered  “implicitly”  varied  by  directly 
replacing  the  velocity  components  in  Equation  68  as  follows: 


C - <n  + |<i|  + «s.>. 

- U!n\\.n  + |fCU  I + «®2>- 

'tn  “ C!  + |C|  <"*,  + -v,). 


and 


(86) 


vi*\ 
m.n - 1 


vi+ 1 
m.n-  1 


+ 


yi+ 1 
m.n - ! 


(wvj 


+ tv  v,). 


provided  „ =/=  0;  i.e.,  provided  an  empirical  tidal  amplitude  is  available  for  the  considered  mesh  cell  Sm  n-  In 
Equations  86,  the  consistency  and  scale  parameters  (u,  u)  and  (v,  v)  are  defined  by 


«1 

= 1, 

Ml 

= 0 

for  ' • Uj*l  <0. 

* m ,n  m.n  * 

M| 

= o, 

Ml 

= A* 
m.n 

otherwise,  but 

Ml 

= 0, 

Ml 

= 0 

* 0; 

(87a) 

Ml 

” 1, 

Ml 

= 0 

f°rAC  • >»• 

Ml 

= 0, 

«2 

= j4* 

m+n.n 

otherwise; 

(87b) 
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, v,  = I.  v,  = 0 


Vj  =0,  v,  = B* 
in.  n 


f°'<xi,  • mi  <o. 


otherwise; 


(87c) 


and 


- v2  = 1,  v2  = 0 


f°'<+in  • vL:\-x  >0- 


vj  = 0,  vj  = B* 


m.n-  1 


v2  = 0,  v,  = 0 


otherwise,  but 


(87d) 


The  continuity  gap  (Eq.  85)  will  be  narrowed  when  the  “control  parameters"  w and  w are  determined  successively  by 

,Af'+1 

" = l 


[Af£, , n forf^O. 


v0  for  f = 0 

with  the  first  “control  limit" 

I tv  I < Ac  i 
and 


(I <Xn  - “fl  /? 


0 


for  f # 0, 
forf  =0 


with  the  second  control  limit 
lit- 1 < Ac2  , 
where  (see  Eq.  68) 


(88a) 


(88b) 


(89a) 


(89b) 


and 


s ■ K i".  K!|  ♦ * v.c*,|C  I * I 

f ■ c'*  i° . te; I * “ > 1 1 1 1 * *•  <■; l c I * >■  c; 


(90) 
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It  is  important  to  note  that  uf  • ^ = 0 and  v(.  • vf.  = 0 for  i = 1,2.  Accordingly,  both  control  limits,  Ac , and  A 2 , 
which  are  at  one's  disposal,  regulate  the  allowed  decrease  or,  respectively,  increase  of  the  velocity  components  in 
Equations  86;  i.e.,  the  implicitly  permitted  corresponding  increase  or  decrease  of  the  local  bottom-friction  coef- 
ficients in  Equations  67.  Since  the  integration  sweeps  across  the  ocean  from  m = to  360  and  n = 4 to  168,  the 
special  choice  of  Uj  = =0  and  v2  = v2  = 0 in  Equations  87a  and  87d  excludes  possible  double  adjustments  of  the 

velocity  components.  Also,  if  u,  * fi",  and/or  v2  * v2,  backward  adjustments  of  the  tidal  elevations  via  the  corre- 
sponding Equation  68  must  be  made.  This  requires  the  replacements 


sm  - n.n 


and 


f. 


/♦I 
m - n,n 


- C1 


Um,\  + > 


C-r CL.  -C|CL.H  +*^)- 

Analogous  substitutions  in  the  forward  directions  of  m and  n follow  automatically  in  the  integration  process. 
The  velocity  replacements  in  Equations  86  may  be  illustrated  by  the  example 

U'*1  > 0,  U>*}  > 0,  f"+1  > 0,  V'*1  , $ 0 

m.n  m*n.n  m.n  ' m.n- 1 < 


(91) 


f 

One  finds  w > 0,  w > 0,  and 


= 0,  £ . #0. 

m-ii.n  ' sm,/i-l 


(92) 


U>*  i - U>*x  (1  + u’/l4  ), 

m,n  m.n  ' m.n’' 


ULXn  - C,,,*1  - ">• 


v>* 1 - vl+l  (i  + tvB4  ), 

m.n  m.n  ' m.n' 


and 


f/+l 

5 m - 


m - p.n 


{>* 1 - C'(//+l  wA* 

- n.n  n m,n  m.n 


(93) 


At  this  point,  it  must  be  mentioned  that  attempts  were  explored  to  lift  the  control  limits  prescribed  by  k ( and 
A2  in  Equations  88b  and  89b  in  an  effort  to  close  the  continuity  gap  completely.  However,  since  the  bottom-friction 
coefficient,  b.  in  A 4 and  B 4 of  Equations  69  is  rather  small  (Eq.  4b),  the  control  limits,  it(  and  k2.  had  to  be  kept 
small  to  achieve  best  results.  Computations  conducted  with  large  control  limits  (excessive  bottom  friction) 
seemed  to  close  the  continuity  gap.  but  the  tidal  and  velocity  fields  in  the  open  oceans  assumed  unrealistically  small 
values  Large  control  limits  k2  (insufficient  bottom  friction)  produced  strong  instabilities  as  anticipated  from  the 
analysis  in  Section  5.G.  To  safely  check  the  possible  instability,  the  second  control  parameter  SP  (Eq’s.  86,  87,  and 
93)  was  defined,  respectively,  in  units  of  5 = ,44  and  ^ = 54  in  contrast  to  u = 1 and  v = 1 , used  for  the  first  control 
parameter  w. 
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(94) 


After  some  trial-and-error  computations,  the  following  control  limits  were  chosen  for  the  M,-tide  model: 

it,  =.03,  k2=  . 06. 

These  moderate  values  reflect  the  well-known  fact  that  the  magnitude  of  bottom  friction  has  a strong  effect  on  the 
motions  considered.  Indeed,  with  some  minor  improvements  of  the  tidal  field,  significant  improvements  of  the 
continuity  gap,  velocity  field,  and  convergence  of  the  integration  were  achieved.  This  procedure  was  applied  to  all 
oceanic  cells  with  known  empirical  tide  data  (Eq.  84),  provided  these  cells  bordered  terrestrial  cells  or  contained 
small  islands  or  other  bottom  irregularities.  No  meaningful  reason  was  seen  to  apply  the  same  bottom-friction 
adjustment  procedure  to  distant  offshore  oceanic  cells  with  available  deep-sea  tide  measurements.  A comparison 
of  computed  tide  data  with  those  empirical  deep-sea  values  is  given  in  Table  4. 

In  order  to  implement  the  second  step  (b)  of  the  hydrodynamical  interpolation  procedure,  the  following 
velocity  replacements  in  oceanic  mesh  cells  bordering  terrestrial  cells  were  defined: 


U'*1  -*  wm ,{/'+' 
m,n  1 m+ii,n 


U'+l  - U'+l 

m+n,n  ^ nt  n 


and 


V>+\  - wv,  K'+»  . 
m,n  1 m,n- 1 


1"+1  , - w V,  V'*1 


m.n-  1 " ' * - m,n 

provided  £ ^ 0 in  Equation  84.  The  parameters  (u  v")are  mutually  consistent  by  definition: 


(95) 


u,  = 1 if = 0,  otherwise  u,  = 0, 
u 2 = 1 if  Uj*}  „ = 0,  otherwise  u2  - 0. 

V,=  1 if  = 0,  otherwise  V,  = 0,  (96) 

and 

V2=  1 if  ^mn - 1 = otherwise  v2  = 0. 

The  remaining  continuity  gap  will  be  further  narrowed  when  the  control  parameter  w is  determined  to  be  in  agree- 
ment with  Equations  68,  85,  88,  89,  and  90  by 

_ -»t-  *?]/?  for?  *0 

10  forf=0  (97a) 
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with  the  third  control  limit 


Iw  I <*j  , 


(97b) 


where 


f = Cn  I"-  Ut:X.n  - “I  O + V.CJC1-, 


V,C3F'+1 
v2  « /n,« 


(98) 


Obviously,  the  substitutions  (Eq's.  95)  specify  consistent  in-  or  out-flows  across  the  mathematical  boundaries 
of  oceanic  coastal  cells  as  illustrated  in  Figure  1 without  explicitly  fixing  the  physical  boundary  line.  Again,  no 
complete  removal  of  the  continuity  gap  was  possible.  The  most  satisfactory  results  for  the  M2  tide  (Part  II)  were 
achieved  by  setting  the  third  control  limit  (Eq.  79b)  at 

*3=0.5.  (99) 

While  the  improvement  of  the  tidal  field  was  again  moderate,  the  remaining  continuity  gaps  and  near-shore  velocity 
distortions  assumed  uniformly  satisfactory  levels.  These  remaining  small  shortcomings  of  the  model  can  easily  be 
attributed  to  the  boundary  inaccuracies  (c),  (d)  and  (e)  listed  above,  but  for  which  no  simple  remedies  were  found. 

It  may  be  emphasized  that  the  rather  significant  change  in  the  near-shore  velocity  field  permitted  by  the  in- 
and  out-flow  specifications  (Eq's.  95)  affected  the  tidal  field  only  in  a minor  fashion.  This  important  phenomenon 
is  in  agreement  with  the  well-known  fact  that  the  pressure  distribution  in  a fluid  motion  is  very  insensitive  to  large 
but  local  velocity  variations.  For  instance,  it  is  perhaps  the  most  important  postulate  in  Prandtl’s  boundary-layer 
theory  (see,  e.g.,  Schlichting.  1968),  and  it  is  the  basis  of  the  hydrostatic-pressure  assumption  invoked  here  in 
Section  3.E  for  the  present  tidal  model. 


G.  STABILITY  ANALYSIS 

A rigorous  stability  analysis  of  the  homogeneous  DOTEs  (Eq’s.  67  and  68)  is,  of  course,  not  possible.  How- 
ever, under  the  assumption  of  constant  coefficients  A,  B,  and  C,  the  simplified  DOTEs  possess  Fourier-type  eigen- 
solutions  (Eq's.  109)  that  permit  a local  stability  analysis  of  the  difference  system.  As  is  well  known  (see,  e.g., 
Richtmyer,  1957),  such  a local  stability  analysis  produces  stability  limits  that  are  usually  sufficient  for  computa- 
tional purposes.  Indeed,  computer  experiments  showed  that  the  stability  limits  so  derived  below  were  scrupulously 
binding  for  the  success  of  the  integration.  The  following  analysis  is  an  expanded  version  of  the  investigation  pre- 
sented by  Zahel  (1970). 

In  detail,  the  following  simplifications  may  be  assumed: 

(a)  b = 0;  i.e.,  no  bottom  friction 

(b)  n = 0;  i.e.,  no  Coriolis  force 

(c)  For  an  arbitrary  but  fixed  mesh  cell  (see  Eq’s.  56,  57,  and  69!,  b), 
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r = rn(«)  = r;|(v)  = 1 /jjsin  6 = locally  constant, 

^ = ^„(w)  = ^„(y)  = 2(I  + 1/P)’ 


w = «n(«)  = wn(v)  =1  + T2; 
and 


(100) 


H = H(u)  = H(y)  = locally  constant, 
m.n  m.n 

H(u)  = H(u)  = H (m)  = 0,(m-v). 
m ,h  m.n  m.n 

It  may  be  mentioned  that  the  assumptions  (a)  and  (b)  have  been  made  in  order  to  display  more  clearly  the  most 
important  stability  characteristics  of  the  DOTEs  that  are  due  to  eddy  dissipation  and  to  the  differencing  parameters, 
k and  k.  It  is  relatively  easy  to  show  that  the  bottom  friction  is  always  stabilizing,  while  the  Coriolii  force  (in  the 
present  differencing  scheme)  is  slightly  destabilizing. 


For  the  following  derivations,  it  is  helpful  to  introduce  some  specified  reference  values  0r  and  Hf  (consistent 
with  Eq’s.  100),  tl^.and  wf,  and 


(101) 


Finally,  the  following  relative  quantities  may  be  introduced: 

£ = li'/i/',.  W = w/wr.  h=H/Hf, 

T=At/£tr.  e=a/ar. 

For  example,  in  this  notation  the  eddy  viscosity,  A (Eq.  6),  assumes  the  form 

A = , 

where  e is  the  “dimensionless  eddy  coefficient." 

The  M ^ -tide  Tables  5-8  (and  Part  II)  are  computed  with  the  reference  values 

0r  = 30°,  //.  = 7 259.84  m, 

so  that 

= 3/4,  ar  = 0.021  919  2 sec’1. 

wr  = 5.  Atr  - 186.309  sec, 

180  oAtrlv  = 1.5°.  Jp  = 360°/ 1 .5°  = 240, 


| (102) 

(103) 


(104) 
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where  J is  the  number  of  time  steps  required  to  integrate  through  one  tidal  period  with  At  = Atf.  Due  to  the 
grading  of  the  grid  system  (Eq.  49),  one  has  almost  everywhere  (n  1,2,3,  and  166,  167,  168) 

1 < £ < j,  1>w>1.  (106) 

Similarly,  due  to  the  cutoff  depth  data  (Eq's.  50,  51 , 71),  one  has  the  limits 

0.001  3 < h < 1 . (107) 

Several  other  reference  values  within  the  stability  limits  have  been  explored.  Particularly  extensive  computations 
were  carried  out  with  Atf  = 248.412  sec  (so  that  J = 180),  but  the  above  reference  values  appeared  to  yield  the 
best  results. 

With  the  simplifications  and  notations  above,  the  coefficients  (Eq’s.  69  and  70)  of  the  homogeneous  DOTEs 
(Eq's.  67  and  68)  become 

A 1 =A2  = Bl  = B2  =0, 

A}  = TB3  = 0GHVAt/L. 

A4  = B4  = 2/ierfco. 

As  =Ab  =BS  =B6  =heT$/ur, 

A1  =A%  =B2  = T2.45, 

A9=Aw  = fl8  =fl9  =B10  =0, 

C1  = VAt/l..  C2  =C3=  At/L. 

The  reduced  DOTEs  (with  constant  coefficients)  yield  the  Fourier-type  eigensolutions 

y>  = Un  jigibjUm  ~ 2/t ) AX  + 7,(2 n - 1)A0), 
m.n  u 1 

yi  = Vnjie>  l7,(2/w  - M)AX  + y (2 n-  1)A0], 
m.n  ° 1 2 

and 

<'  = t„</V  l7,(2»»  - /'(AX  + 7 (2n  - 1 ) AO ] 

50  1 1 

with  an  arbitrary  wave  vector  (7p  >2)  a,'d  50,116  nonzero  amplitude  vectors  (t/Q,  FQ,  f0),  provided  the  eingenvalue 
d satisfies  the  cubic  characteristic  equation 

Mu  - dAi0)  0 

I ' 

0 (A  22  — dA  io ) 

(1  — k + Kd)An  (1  - k + Hd)  A 32 
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A 13 

A 23  =0,  (110) 

(1  -d) 


where  (after  some  algebra) 


L 


1 


/4io  = 1 + kA*  = 1 + 2k her ^ to  , 

A 1 1 = A 22  = A io  — 2A 4 s2 , 
and 

^13^31  + -4  32-4  23  = - 4(3hwr2SJ 

with 

, T2  sin27.A0  + sin2>2^0 
0 < s2  < 4 - < 1. 

r2  + 1 

The  cubic  characteristic  in  Equation  1 10  yields  the  three  eigen1  aiues  (dQ,  d = d,,  d 2) 


’ 

(111a) 

(111b) 


and 

dQ  = 1 - Aher^icjs2/A  j , 

i r - -.fi1/2 

(112a) 

dA  10  = A 10 

where 

- 2Ar«sM0,  ± 2 its  huAfiA  ,Q  - /icosMq,)  , 

(112b) 

The  DOTEs  will  be  stable,  provided 

-4  oi  =eif  + 0KT. 

(112c) 

Idj<  1 for  It  = 0, 1,2. 

(113) 

Under  the  strict  inequality  of  Equation  1 1 3,  the  three  eigenvalues  dQ,  d , , and  d2  define  three  decaying  eigenwaves 
represented  by  Equations  109.  Since  dQ  is  real,  the  corresponding  eigenwave  is  a standing  wave  with  no  phase  shift 
if  d0  > 0 (see  Eq.  1 19).  The  other  two  eigenvalues  d,  and  d2  define  a pair  of  eigenwaves  progressing  in  opposite 
directions  with  the  same  decay  and  dispersion  rates,  provided 

PAw>has2A20V  (114) 

This  condition  holds  true  for  all  0 < s <1  (Eq.  111b);  i.e.,  for  all  wave  vectors  (y, , y2)  if  and  only  if 

(Hjo^hC^.  (115) 

If  this  condition  fails,  then  there  exist  some  short  waves  with  large  wave  numbers,  y,  and  y2,  which  become  stand- 
ing waves  of  different  decay  rates.  However,  all  sufficiently  long  waves  remain  dispersively  progressing  and  decaying 
at  the  same  rates. 

It  seems  physically  plausible  to  treat  all  long  and  short  waves  equally  and.  hence,  to  impose  the  conditions  of 
Equations  1 13  and  1 15  on  the  free  parameters,  e,  t,  k,  and  k.  Using  Equations  1 12, 1 13,  and  115,  one  finds 

Id,  I2  = Idjl2  * 1 - 4hrC3s2  [ej’-(l  -k)0t]  />4,0  (116) 

with 

Id,  I2  - ld2l2=>d0forK  = l.  (117) 
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For  x < 1.  the  stability  condition  requires 


I 


I 


' 


-k)0t;  (118) 

i.e.,  a minimum  of  eddy  viscosity  is  necessary  for  stability.  However,  for  k = 1 , no  minimum  eddy  viscosity  is  re- 
quired, wliich  explains  the  choice  made  here  (Eq.  72)  and  by  Hansen  (1966),  Zahel  (1970,  1973,  1975),  and  Estes 
(1975,  1977). 


For  the  chosen  value  k = 1 , the  stability  condition  is  satisfied  for  all  s and  all  eigenvalues  <i0 . d, , and  d2  when 


4herifu) 

(119) 

i.e.,  when 

r = < 1/2(2  - x)eliiptu, 

(120) 

provided  (Eq.  115  in  explicit  form)  also 

> 02t2  + 20er^  (1  - k)  + e 2<J/2. 

(121) 

hui 

The  obviously  increased  stability  limits  unposed  by  both  Equations  120  and  121  explain  the  choice  of  x = 1 made 
here  for  the  present  tide  model  (Eq.  72)  in  deviation  from  the  value  k = 0 used  by  Hansen,  Zahel,  and  Estes.  It  may 
also  be  recalled  that  as  an  important  by-product,  k = 1 facilitated  the  simple  hydrodynamical  interpolation  of 
empirical  tidal  data  into  the  model  described  in  Section  5.F. 

With  (3  = 0.90  (Eq's.  17)  and  the  possible  extreme  values  of  h = 1 and.  simultaneously,  if  = tc  = 1 (Eq's.  106 
and  107),  one  finds  from  Equations  120  and  121  for  k = 1 and  t = 1 the  allowable  range  for  the  dimensionless 
eddy  coefficient,  e . 

0 < e < 0.3  (At  = Atr).  ((22) 

The  same  range  holds  also  for  the  southern  three  colatitude  lines  n = 166,  167,  168,  which  violate  the  condition  of 
Equation  49,  but  for  which  the  relative  depth,  h.  falls  sufficiently  below  unity.  The  upper  limit  on  e could  be 
raised  somewhat  by  considering  the  simultaneous  values  of  u>,  and  h on  each  colatitude  r = 4,  5 168  sepa- 

rately. In  order  to  obtain  the  best  possible  tidal  field,  extensive  trial-and-error  computations  led  to  the  choice 

Ar  = Atr  = 186.309  sec  and  e = 0.075.  (123) 

This  final  choice  completes  the  detailed  parameter  specifications  of  the  present  tide  model. 

At  this  point,  it  may  be  noticed  that  the  stability  requirement  for  the  DOTEs  restricts  the  possible  amount 
of  eddy  dissipation.  As  is  physically  plausible,  the  finite-differencing  parameters,  k and  IT,  the  mesh  si u,  L,  the 
depth,  H,  the  time  step,  A r,  and  the  dimensionless  eddy  coefficient,  e,  are  intimately  related  to  each  other.  Trial- 
and-error  computations  are  needed  to  select  those  parameters  for  best  results.  It  is  particularly  important  to  observe 
that  (especially  for  if  = 1)  the  rate  of  decay  (Eq’s.  112a,  116,  and  1 17)  of  all  eigenwaves  depends  directly  on  the 
product  he.  Accordingly,  for  fixed  e,  waves  in  deep  ( h =»  1)  ocean  basins  decay  faster  than  in  shallow  (h  < 1)  regions 
if  bottom  friction  is  negligible. 
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It  is  obviously  this  physically  realistic  phenomenon  that  led  to  the  introduction  of  the  novel  depth -dependent 
eddy  viscosity,  A,  defined  by  Equations  6 or  103.  For  a depth-independent  eddy  viscosity,  one  has  he  = constant, 
in  which  case  waves  would  decay  at  the  same  rate  in  deep  or  shallow  (see  Sect.  6.A)  oceans,  even  though  no  bottom 
friction  is  present.  Following  Hansen  (1966),  Zahel  (1970,  1973,  1975).  and  Estes  (1975,  1977),  the  present  tide 
model  also  used  at  first  a constant  eddy  viscosity  with  rather  disappointing  results  caused  by  the  strongly  varying 
bathymetry. 

It  is  interesting  to  note  that  for  the  limiting  case  of  Equation  1 1 5 , i.e.,  for 

PAlQ=hZA20l,  (124) 


Equation  1 1 2 assumes  the  simple  form 


dAm  =Am  - 20ts2  ± 2/0rs(l  -s2)* 


Hence,  for  the  northeast  waves  under  45°  with  wave  numbers  (Eq.  1 1 lb) 


2 2 

= 72  = 7,  s -sin  7A0, 


and  (k  = k = 1) 


d = e*  + ^e±2-Ae  =d  + de±2 
e\l/  +0r 

Hence,  in  this  case,  the  eigenvalues and  c/2  lie  on  the  circle 


\d-dc\  = dr,  dc+dr=  1 


as  illustrated  in  Figure  5. 


Figure  5.  Illustration  of  Eigenvalues 
dQ,  d j , and  d2  = tT, , in  Circle 
I d-d  \<d 

c r 
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6.  DISCUSSION  OF  THE  TIDE  MODEL 


A.  DISCRETE  VERSUS  CONTINUOUS  OCEAN-TIDE  EQUATIONS 

In  Section  2(g),  it  was  contended  that  discrete  ocean-tide  equations  (DOTEs,  Eq’s.  67  and  68)  reflect  the 
physical  reality  of  ocean  tidal  currents  more  perfectly  and  perceptibly  than  the  corresponding  continuous  equations 
(COTEs,  Eq’s.  29  and  30).  Although  the  latter  follow  from  the  former  by  a “formal”  limit  process,  it  is  neither 
technically  feasible  nor  theoretically  desirable  to  seek  convergence  of  the  discrete  solution  to  the  continuous  inte- 
gral. In  fact,  in  fluid-flow  problems  of  global  dimensions,  it  is  not  possible  to  approximate  the  conditions  of  the 
continuous  case  to  any  reasonable  degree.  Even  with  future  computer  technology,  it  will  not  be  meaningful  to 
refine  significantly,  for  instance,  the  1°  by  1°  grid  system  defined  in  Sect.  5.A,  which,  with  its  100-km  mesh  size,  is 
far  from  being  anywhere  near  a continuous  description.  Any  attempt  to  refine  the  grid  system  would  have  to  be 
matched  by  an  improved  bathymetry  (Sect.  5.B),  which  requires  worldwide  in  situ  measurements. 

From  the  theoretical  point  of  view,  it  is  equally  superfluous  to  seek  an  e-approximation  of  the  continuous 
situation,  which  in  fact  is  only  vaguely  defined  (Sect.  3. A).  In  laminar  viscous-flow  theory,  such  notions  as  particle 
(point)  velocity  and  pressure,  as  well  as  the  Navier-Stokes  equations,  are  all  derived  from  physically  sound  discrete 
(finite-difference)  definitions  by  a formal  limit  procedure  ; i.e.,  by  simply  assuming  the  existence  of  the  limit  values 
(see,  e.g.,  Schlichting,  1968;  and  Whitaker,  1968).  While  this  assumption  is  well  justified  in  most  laminar-flow 
problems,  it  is  well  known  today  (see.  e.g.,  Ladyzhenskaya,  1969)  that,  in  general,  even  laminar  motions  must  be 
sought  in  the  class  of  generalized  (distribution)  functions.  Hence,  velocities,  pressure,  and  their  derivatives  do  not 
exist  in  the  ordinary  sense  (pointwise);  only  their  “functionals”  (effects  such  as  mass  fluxes,  forces,  momenta)  are 
physically  defined. 

The  ambiguity  of  continuous-flow  models  becomes  much  more  apparent  in  the  critical  laminar  regime.  Experi- 
ments (e.g.,  Busse  and  Whitehead,  1971)  and  theory  (e.g.,  Schwiderski,  1972)  have  clearly  established  that,  when 
given  characteristic  flow  parameters  (dimensions,  velocities,  etc.)  exceed  certain  critical  values,  the  corresponding, 
uniquely  existent  laminar  motions  become  unstable  and  bifurcate  into  laminar  flows  of  infinitely  many  different 
shapes.  The  classical  laminar  boundary  and  initial  conditions  (see,  e.g..  Sect’s.  3.B  and  5.E)are  no  longer  sufficient 
to  specify  a unique  motion.  The  situation  seems  to  be  governed  by  hysteresis  and  pure  chance  rather  than  rigorous 
physical  selection  principles. 

Critical  laminar  motions  become  still  more  microscopically  undefinable  when  the  corresponding  characteristic 
flow  parameters  (as  the  global  dimensions  of  ocean  currents)  exceed  further  supercritical  points  and  the  motions  go 
turbulent.  The  statistical  approach  underlying  the  “time-averaging”  process  to  derive  the  so-called  Navier-Stokes 
equations  of  mean  turbulent  flow  (Eq’s.  1 and  2)  is  entirely  formal  and  vague  (see.  e.g.,  Schlichting,  1968).  For 
example,  what  velocity  (particle,  point,  etc.)  is  averaged  over  what  time  interval?  In  this  respect,  the  present  periodic 
tidal  motions  are  clearly  most  illuminating  of  the  problems  at  hand.  If  one  averages  (as  usually  meaningful)  over 
a “sufficiently"  long  time  span  (say,  longer  than  the  tidal  periods),  then  the  averaged  velocity  should  approach  zero, 
which  is  obviously  not  of  interest. 

Evidently,  turbulent  particle  velocities  manifest  themselves  statistically  through  their  integrated  (macroscopic) 
physical  effects,  such  as  mass  fluxes.  Hence,  the  proper  mathematical  representation  of  turbulent  motions  should  be 
sought  in  the  class  of  generalized  functions.  Since  the  product  of  generalized  functions  has  no  mathematical  meaning 
(see,  e.g.,  Shilov,  1968),  it  appears  understandable  that  there  is  no  way  to  define  the  Reynolds  stress  tensor  of 
turbulent  motion  (Sect.  3.C)  by  a meaningful  ordinary  or  generalized  function,  because  it  contains  quadratic  prod- 
ucts of  the  so-called  fluctuating  velocity  residuals  (see,  e.g.,  Schlichting,  1968).  However,  its  energy-dissipating 
(stress-like)  effect  is  physically  quite  apparent  and  must  be  modeled  in  some  macroscopic  sense. 
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To  avoid  all  conceptual  difficulties  of  microscopic  turbulent  motions,  it  seems  natural  to  fall  back  to  the 
discrete  (macroscopic)  description  of  laminar  flows  that  leads  formally  to  the  Navier-Stokes  equations  (see,  e.g., 
Schlichting,  1968;  and  Whitaker,  1968).  In  fact,  by  proper  “generalized”  interpretation,  virtually  every  notion  used 
in  the  laminar  regime  retains  its  physical  meaning  in  the  discrete  turbulent  domain.  For  example,  the  “mean” 
x -velocity,  u,  of  a "flow  parcel”  contained  in  a rectangular  test  (grid)  cell  of  mesh  lengths  Ax,  Ay,  and  A z (Figure 
6a)  at  some  time,  /,  of  a time  interval.  At,  is  defined  as  the  mass  flux  (AMX)  crossing,  say,  the  central  (x  = Xj) 
surface  element  (Ay,  Ac)  of  the  cell  during  the  time  span  Ar  divided  by  the  fluid  density, p,  the  area,  Ay  • Az,  and 
the  time.  At,  so  that 


u = AMX  IpAyAxAt. 


(129) 


AM*  = upAyAzAt 


AX/2  Xf  AX/Z 


Figure  6a.  Illustration  of  Average  Velocity; 

AMX  - Mass  Flux  in  x-Direction, 
u = Mean  Velocity  in  x-Direction, 
and  p = Fluid  Density 

This  (generalized)  definition  of  the  mean-flow  velocity,  u,  is  uniformly  valid  for  laminar  and  turbulent  mo- 
tions, since  A Mx  is  a physically  realistic  (measurable)  quantity  for  any  p > 0,  Ax  > 0,  Ay  > 0,  Az  > 0,  and  At  > 0. 
If  u exists  in  the  limit  Ax  -*■  0,  Ay  -*■  0,  Az  -*•  0,  and  Ar  -*  0,  then  u becomes  the  ordinary  particle  (or  point)  velocity. 
This  assumption  is  well-justified  in  most  (see  above)  laminar  motions,  where  neignboring  particles  flow  on  smooth 
neighboring  pathlines.  In  contrast  to  laminar  flow,  the  limit  assumption  is  certainly  not  justified  in  turbulent  flow, 
where  the  fluid  particles  move  in  an  indeterministically  “turbulent,”  “fluctuating,”  or  “eddying”  way.  In  the  sense  of 
generalized  functions,  u could  be  defined  in  the  limit  by  prescribing  its  mass  flux  (A Mx)  for  every  nonvanishing  p, 
Ax,  Ay,  Az,  and  At. 

Using  similarly  generalized  definitions  of  all  mean  velocity  components  and  pressure,  and  retaining  their  finite 
differences  in  place  of  (generally  non-existentjordinary  derivatives,  one  arrives  as  usually  at  the  discrete  Navier-Stokes 
equations  (analog  of  Eq’s.  1 and  2)  without  requiring  the  existence  of  any  limit  values  (see,  e.g.,  Schlichting,  1968; 
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and  Whitaker,  1968).  They  express  simply  the  physical  laws  of  conservation  of  momentum  and  mass  for  every  test 
(mesh)  cell.  This  is  particularly  tangible  for  the  DOTE  (Eq.  68),  which  conserves  mass  by  balancing  the  excess  of 
mass  flux  into  a mesh  cell  with  a corresponding  increment  in  tidal  height. 

In  the  discrete  form,  the  Navier-Stokes  equations  look  formally  the  same  for  laminar  and  turbulent  motions. 
Yet,  when  in  the  laminar  case  velocity  and  pressure  exist  pointwise  in  the  limit,  one  has  a unique  microscopic 
(continuous)  description  of  the  flow  subject  only  to  specified  boundary  and  initial  conditions.  Hence,  the  quality 
of  a solution  of  the  discrete  model  can  be  measured  against  the  so-called  exact  integral  of  the  continuous  equations. 
In  sharp  contrast  to  the  laminar  case,  no  unique  microscopic  description  of  turbulent  flow  exists  against  which  he 
quality  of  the  discrete  model  could  be  measured.  Therefore,  the  numerical  analyst  must  seek  an  optimum  size  of 
the  mesh  cell  (Ax,  Ay,  Az)  and  the  time  span  (Ar)  in  order  to  model  the  undetermined  microscopic  turbulent 
motion  so  that  their  macroscopic  effects  match  the  expected  or  observed  features. 

The  strong  dependence  of  a discrete  turbulent-flow  model  on  the  size  of  the  mesh  cell  (Ax,  Ay,  Az)  and  the 
time  span  (Ar)  can  be  assessed,  for  instance,  from  the  definition  of  the  average  velocity,  u,  by  Equation  12°.  With 
increasing  mesh  area  (Ay  • Az),  more  and  more  as  well  as  larger  and  larger  fluctuating  or  eddying  motions  are 
filtered  out  and  remain  unaccounted  for  in  the  average  value  of  u.  Hence,  the  maximum  mesh  lengths  Ax,  Ay,  and 
Az  must  e sufficiently  smaller  than  the  smallest  wave  length  one  wishes  to  resolve.  On  the  other  side,  if  the  area 
(Ay  • Az)  in  Equation  129  is  chosen  smaller  and  smaller,  then  u becomes  more  and  more  undetermined  (fluctu- 
ating). 

Similar  arguments  determine  an  optimum  time  step,  Ar.  In  the  present  discrete  tide  model,  the  cell  size  was 
reasonably  limited  by  the  available  bathymetric  tables.  The  time  step  Af  (Eq.  1 23)  was  determined  by  trial-and-error 
compulations  so  that  60  time  points  represent  one-quarter  of  the  M^-tide  period. 

Another  significant  distinction  between  the  discrete  Navier-Stokes  equations  of  laminar  and  turbulent  flow 
becomes  apparent  in  the  average-stress  tensor.  As  is  well  known,  the  turbulent  fluctuations  neglected  in  the  mean 
velocity  and  pressure  manifest  themselves  as  stress-like  (energy -dissipating)  forces  that  affect  the  mean  motion. 
Unfortunately,  no  exact  and  unique  constitutive  equation  is  known  today  that  relates  those  turbulent  Reynolds 
stresses  to  the  mean  rate  of  strain  (deformation  of  the  flow  parcel)  determined  by  the  average  velocity  (see,  e.g., 
Schlichting.  1968;  and  Whitaker,  1968).  In  the  Boussinesq  (1877)  substitution  used  in  the  present  tide  model,  the 
mean  turbulent-stress  tensor  is  directly  related  to  the  average  rate  of  strain  (analog  of  Eq's.  5).  Hence,  the  macro- 
scopic stress  eftects  of  the  turbulent  fluctuations  on  the  mean  velocity  are  assumed  to  follow  a similar  simple  law 
as  the  viscous  laminar  stresses  in  Newtonian  fluids.  Only  the  coefficient  of  viscosity  is  replaced  by  the  so-called 
edd-'  viscosity  which  remains  to  be  modeled  to  account  for  the  otherwise  neglected  eddying  motions  in  some  best 
sense. 


In  the  absence  of  better  approximations,  it  seems  idle  to  argue  endlessly  about  the  physical  justification  of 
the  Boussinesq  substitution;  the  fact  remains  that  it  represents  the  simplest  possible  constitutive  equation,  including 
zero  used  by  some  researchers.  Moreover,  it  provides  for  considerable  flexibility  to  model  the  microscopically 
undetermined  but  maeroscopically  apparent  eddy  dissipation  by  choosing  suitable  velocity-dependent  or  velocity- 
independent  eddy  viscosities  cither  uniformly  or  separately  for  all  three  stress  directions.  Evidently,  if  the  velocity 
field  were  known,  a priori,  then  one  could  determine  exact  eddy  viscosities,  a posteriori,  in  many  ways.  In  this 
connection,  it  is  of  interest  to  know  that  the  mean  flow  is  quite  insensitive  to  fairly  large  variations  (say  25  percent) 
of  the  eddy  viscosity.  This  observation  by  Munk  and  Palm^n  (1951)  was  confirmed  for  oceanic  tidal  motions  by 
the  author's  extensive  computer  experiments.  It  is  probably  related  to  the  well-known  fact  that  potential  motions 
satisfy  the  complete  Navier-Stokes  equations  of  laminar  flow  with  any  constant  viscosity.  Above  all.  as  with  any 
other  physical  law,  the  Boussinesq  substitution  has  successfully  passed  its  crucial  test  in  many  practical  applications 
in  hydrodynamics,  oceanography,  and  meteorology.  The  present  ocean-tide  model  is  no  exception  (see  Sect.  6.B). 
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In  order  to  illustrate  the  Boussinesq  substitution  in  the  discrete  case,  one  may  consider,  for  example,  the 
average  normal  stress  r**  produced  by  the  filtered  out  (see  the  remarks  to  Eq.  129)  fluctuating  motions  on  the 
surface  (Ay,  A z)  at  x = x(  shown  in  Figure  6b.  Following  Boussinesq.  one  has  (see  the  analogous  Eq.  5c) 

rf  = 2Ap^--^,  (130a) 

1 2 Ax 

where  uQ  and  «2  are  the  corresponding  mean  velocities  at  x = xQ  = x(  Ax  and  x = x2  = Xj  + Ax.  Hence,  the 
turbulent  stress  r™  grows  linearly  with  the  rate  of  change  of  average  velocity.  Analogous  to  the  corresponding 
laminar  stress,  this  linear  law  appears  physically  acceptable,  since  the  expected  mean  tidal  velocities  are  small.  One 
concludes  from  Equation  130a  that  a large  change  of  mean  flow  produces  a large  turbulent  stress  which  plausibly 
must  be  due  to  strong  fluctuating  motions. 


Figure  6b.  Illustration  of  Mean  Normal  Stress: 

Mq,  u2  = Average  x-Velocities  at  x = xQ.  x2. 
and  rj  ' = Average  Normal  Stress 

If  one  assumes  a constant  eddy  viscosity.  A.  in  Equation  130a.  then  the  analogy  between  turbulent  and 
laminar  stress  becomes  complete.  However,  as  was  pointed  out  above,  the  strength  of  the  fluctuating  motions 
under  consideration  depends  on  the  size  of  the  surface  area  (Av  • Ac),  which  justifies  the  assumption 

A = aAyAz  (130b) 

in  Equation  130a.  where  a may  be  held  constant  or  used  for  further  modeling.  Similar  arguments  can  be  developed 
for  all  six  turbulent  stress  components  (see  Eq's.  5)  with  eddy  viscosities  equivalent  to  Equation  1 30b. 

The  novel  eddy-viscosity  law  expressed  by  Equation  130b  is  obviously  equivalent  to  the  eddy  viscosity  intro- 
duced in  the  present  tide  model  by  Equation  6a.  It  also  explains  Equation  4a,  which  specifies  the  bottom-friction 
coefficient,  B,  of  the  discrete  tide  model.  The  need  for  a mesh-area-dependent  eddy  viscosity  became  apparent  when 
initial  tidal  computations  with  a constant  value  failed  to  yield  realistic  results  (see  Sect's.  2(a)  and  5.G).  It  must  be 
emphasized  that  the  microscopically  indeterministic  nature  of  turbulent  motions  is  not  completely  removed  in  the 
discrete  flow  model.  Its  specific  macroscopic  effects  on  the  mean  flow  are  apparent  in  the  required  optimum  choice 
of  the  grid  system,  time  step,  differencing  parameters,  and,  most  of  all,  the  eddy-viscosity  law  and  scaling  coeffi- 
cient. No  information  on  the  fine  structure  of  the  turbulence  can  be  expected  from  such  a model. 
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B.  QUALITY  OF  THE  OCEAN-TIDE  MODEL 


Because  of  the  absence  of  an  exact  (continuous)  ocean-tide  model,  the  degree  of  reality  achieved  by  any 
approximate  (discrete)  model  must  be  measured  against  empirically  known  features  of  ocean  tides.  In  contrast  to 
other  turbulent-flow  problems,  for  instance  in  general  ocean  currents,  a large  number  of  tidal  observations  (Sect. 
5.C)  around  the  world  are  available  for  comparison.  Since  the  present  tide  model  incorporates  essentially  all  known 
empirical  data  by  hydrodynamical  interpolation  (Sect.  5.F),  no  direct  comparison  of  observed  and  computed  data  is 
feasible. 

Nevertheless,  a rather  vivid  and  comprehensive  appraisal  of  the  reality  of  the  present  tide  model  is  possible  by 
inspecting  the  quality  of  hydrodynamical  interpolation;  i.e.,  by  evaluating  the  "smoothness”  with  which  the  com- 
puted tide  "accepts  or  rejects”  (Figure  7)  the  empirical  tidal  data.  In  fact,  the  smoothness  characteristics  of  the 
novel  hydrodynamical  interpolation  technique  are  distinctly  different  from  those  of  other  direct  interpolation  pro- 
cedures using  power  or  trigonometric  polynomials.  In  the  latter  case,  smoothness  of  the  interpolation  can  be  carried 
up  to  any  desired  degree  by  simple  intend.  The  adjustment  of  hydrodynamical  parameters  (Sect.  5.F)  in  the  former 
method  does  not  imply  any  smoothness  of  the  interpolation,  unless  both  the  empirical  input  data  and  the  hydro- 
dynamical tide  model  are  compatible  with  each  other.  As  is  well  known  (see  Sect’s.  S B and  5.C).  local  tidal  distor- 
tions, caused  by  an  isolated  roughness  (seamount  or  small  island)  in  the  bottom  relief,  affect  the  surrounding  ocean 
tide  very  little.  The  major  level'of  ocean  tides  is  shaped  by  continental  shorelines  and  large  (in  area  and/or  length) 
islands  and  ridges.  In  contrast  to  ordinary  polynomial  interpolations,  an  important  feature  of  the  new  hydrody- 
namical interpolation  method  is  that  it  preserves  those  significant  properties  of  ocean  tidal  currents  without  any 
essential  alterations. 


Figure  7.  Illustration  of  "Rejected”  (Protruding)  Empirical  Tidal  Amplitude 
or  Phase  Value:  x = Rejected  Input  Tide  Values,  and 
® = Computed  Tide  Values 

Extensive  computer  experiments  were  conducted  to  test  the  important  smoothness  characteristics  ol  the 
hydrodynamical  interpolation  procedure.  Faulty  input  data  were  deliberately  inserted  and  quickly  recognized  as 
rejected  by  (protruding  out  of;  Figure  7)  the  computed  surrounding  tide.  Indeed,  the  first  computations,  which 
included  empirical  tidal  data,  revealed  immediately  several  input  errors  in  the  data.  Vice  versa,  smoothly  accepted 
empirical  tidal  data  were  randomly  deleted  to  test  their  backlash  reaction  on  the  computed  tide.  As  anticipated,  no 
significant  modifications  were  detected.  Consequently,  the  hydrodynamical  interpolation  technique  permits  a check 
of  the  reality  of  both  the  tide  model  and  the  empirical  tidal  input  data.  If  an  input  value  is  rejected  by  the  computed 
tide,  then  one  or  the  other  or  both  are  defective.  Fortunately,  only  very  few  discrepancies  between  the  different 
sources  of  observed  Mj-tide  data  (see  Sect.  5.C)  have  been  discovered  that  way. 


The  present  discrete  tide  model  has  been  applied  to  compute  the  global  M-,  ocean  tide.  A complete  discussion 
and  tabulation  of  all  amplitudes  and  phases  will  be  presented  in  Part  II  of  this  report.  In  order  to  display  the  quality 
of  the  tidal  model,  the  computed  amplitudes  (in  cm)  and  phases  (in  degrees)  along  with  their  adjacent  empirical 
values  have  been  tabulated  in  “30°  by  50°  map  form"  for  four  typical  ocean  areas  (Tables  5-8).  All  empirically 
supported  input  data  along  continental  shores  and  at  island  stations  are  underlined  in  the  tables.  All  near-shore 
deep-sea  measurements  included  in  the  model  are  labeled  by  subbrackets.  As  was  explained  in  Section  5.F,  all 
distant  offshore  deep-sea  measurements  are  not  included  in  the  tide  model.  However,  their  approximate  locations 
are  marked  by  wavy  underlines  and  their  corresponding  observed  data  are  listed  in  Table  4.  Land  points  are  left 
blank. 


In  the  evaluation  of  the  tidal  accuracy,  one  must  remember  that  the  ocean  tide  at  any  fixed  location  is  deter- 
mined by  two  harmonic  constants.  If  (£0,  60)  and  (£,  6)  denote  the  respective  local  amplitudes  and  phases  of  the 
“true"  and  "computed"  tides 


f0  = £0cos(or-50),  f = S cos  (or  - 6), 

(131) 

then  their  time-dependent  error  is 

F=f0  = £cos(or  5) 

(132a) 

with  the  standard  deviation 

rms(f)=7V/?T- 

(132b) 

where 

?2=So  2{0«cos(So  6)  + {2 

(132c) 

and 

S0  sin  60  ( sin  5 

Ian  5 = z r 

S0cos5o  «cos6 

(I32d) 

Some  maximum  erros  arc- 

W«0  + * for  5o-fi  = l80°’ 

(133a) 

f w = l;M  - £o  “ ? f°r  60  - 8 = 0 . 

(133b) 

7„  = sinT(6o  “ 5)  for  * = $<r 

( 1 34a) 

and 

TM=TM=i  for  £ = s0  and  80  -6  =60° 

(134b) 

Equation  134b  expresses  the  important  fact  that  a 60°  phase  error  results  in  an  amplitude  error  equal  to  the  tidal 
amplitude  and.  hence,  renders  the  computed  tidal  prediction  completely  useless.  Of  course,  in  regions  of  sufficiently 
small  amplitudes,  any  phase  error  is  acceptable. 

Tables  5a  and  5b  depict  the  tidal  amplitudes  and  phases,  respectively,  of  the  northwestern  Atlantic  Ocean 
including  the  eastern  Caribbean  Sea.  As  can  be  verified  by  earlier  tide  models,  this  entire  area  was  very  difficult  to 
model,  because  its  rough  bottom  topography  has  a strong  effect  on  the  tidal  currents  that  sweep  over  or  across 
various  barriers  with  rapidly  changing  water  levels.  There  is  the  broad  and  shallow  continental  shelf  along  the  whole 
North  American  shoreline  with  Cape  llatteras.  Long  Island.  Cape  Cod.  Nova  Scotia,  and  Newfoundland  all  protrud- 
ing into  the  ocean  basin.  Furthemiorc,  there  are  the  Grand  Banks,  the  Bahama  Banks,  and  the  long  and  narrow 
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Caribbean  Ridge.  Obviously,  all  of  the  corresponding  local  tidal  features  could  not  be  realistically  captured  by  the 
tide  model  without  a proper  representation  of  the  bathymetry  (Sect.  5.B)  and  without  the  hydrodynamical  inter- 
polation (Sect.  5.F)of  the  locally  collected  tidal  observations. 


Now,  if  one  scans  the  tidal  amplitudes  and  phases  (Tables  5a  and  5b)  from  the  north  to  the  south,  one  gathers 
the  vivid  impression  that  the  whole  computed  ocean  tide  is  completely  locked  into  the  array  of  empirical  (under- 
lined) tidal  data  everywhere  along  the  continental  coast  and  along  the  many  aligned  islands  separating  the  Atlantic 
Ocean  from  the  Gulf  of  Mexico  and  the  Caribbean  Sea.  It  is  particularly  impressive  to  see  the  observed  tide  data  at 
the  offshore  islands  (Sable-Si,  Barbados-BB.  and  even  as  far  as  Bermuda-Bl)  and  at  the  included  near-shore  (sub- 
brackets),  deep-sea  stations  all  realistically  well-accepted  by  the  computed  surrounding  tide.  Moreover,  one  finds 
the  excluded  offshore  deep-sea  measurements  (locations  marked  by  wavy  underlines)  in  the  Atlantic  and  Caribbean 
Sea  fully  verified  by  the  independent  tide  model. 

As  can  be  seen  in  the  special  listing  of  Table  4,  the  measured  and  computed  amplitudes  and  phases  at  the 
Atlantic  stations  agree  within  2 cm  and  6°.  respectively.  The  remaining  discrepancy  is  probably  within  the  experi- 
mental error  due  to  short  observation  times  and  the  use  of  the  distant  reference  station  Bermuda  (Zettler  et  al., 
1975),  which  exhibits  even  larger  gaps  between  the  various  tidal  observations  listed  in  Table  3. 

Attention  may  be  drawn  to  the  existence  of  considerable  slopes  between  the  empirical  boundary  data  and  the 
computed  ocean-tide  values  in  the  high-amplitude  ranges  from  Nova  Scotia  to  Cape  Cod  and  from  Cape  Hatteras 
to  Florida’s  coast.  Yet.  these  rapid  tidal  variations  can  be  considered  as  realistic  because  throughout  the  same  sec- 
tions the  empirical  data,  amongst  themselves,  display  exactly  the  same  roughness.  This  only  substantiates  clearly 
the  fundamental  difference  between  polynomial  and  hydrodynamical  interpolation  techniques  pointed  out  above. 

In  Part  II  of  this  report,  the  same  tidal  roughness  will  be  recognized  in  several  similar  coastal  places  around 
the  world.  From  this  typical  phenomenon,  one  can  draw  the  fortunate  conclusion  that,  while  some  empirical  data 
may  be  lacking  high  accuracy  (see  Table  2 and  the  British  Admiralty  Tide  Tables,  1977),  the  computed  adjacent 
ocean  tide  may  retain  its  high  quality. 

In  order  to  gain  a deeper  insight  intc  the  detailed  tidal  phenomena  from  the  enclosed  table  charts  (say.  Tables 
5a  and  5b).  it  is  helpful  to  recall  the  physical  meaning  of  the  tabulated  tidal  constants.  The  local  tidal  amplitude,  £, 
is  defined  as  half  the  tidal  "range,"  which  measu.-s  the  total  variation  of  the  water  level  from  high  to  low.  Lines  of 
constant  amplitudes  are  called  "corange  lines.”  As  can  be  seen  from  Equations  12  through  15,  with  X = 0,  the 
local  phase.  6,  specifies  the  tidal  cresting  time  (in  degrees)  after  the  moon's  (or  sun’s)  passage  over  the  Greenwich 
meridian  (X  = 0).  For  the  present  M2  tide  (see  Table  1),  one  has  the  following  time  conversions: 

360°  = 12.421  hours  (period) 

30°  = 1.035  hours 
1°  = 2.070  minutes 

Lines  of  constant  phases  (simultaneous  cresting  times)  are  called  “cotidal  lines.”  In  particular,  at  the  0°  = 360° 
cotidal  lines,  which  are  conspicuously  visible  in  the  phase  charts  (Tables  5b  to  8b),  the  tide  crests  simultaneously 
with  the  moon's  passage  over  the  Greenwich  meridian.  The  tidal  crest  advances  with  time  normal  to  the  cotidal 
lines  toward  larger  phases.  A point  of  zero  amplitude  (J  = 0)  around  which  the  tidal  crest  rotates  from  0°  to  360° 
is  called  an  “amphidromic  point;”  it  is  marked  in  the  tables  by  a circled  star®. 

In  the  area  of  Tables  5a  and  5b,  a major  amphidromic  point  is  visible  in  the  Caribbean  Sea  southeast  of  the 
island  of  Puerto  Pico  (PRI)  near  the  marked  deep-sea  gauge  station.  The  loosely  connected  Caribbean  and  Atlantic 
tides  rotate  counterclockwise  around  this  point  with  the  0°  = 360°  cotidal  line  running  northeastward.  As  a result 
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Table  5b.  M2  Tidal  Phases,  6 (°),  of  the  Northwestern  Atlantic  Ocean 
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of  this  rotation,  the  whole  Caribbean  Sea  appears  to  be  trapped  and  unable  to  develop  any  significant  M2  tide.  In 
agreement  with  the  observations,  the  M2  tidal  crest  sweeps  across  the  Caribbean  Sea  essentially  from  north  to 
south  with  very  little  variation  in  water  level. 

If  one  follows  the  tidal  crest  around  the  amphidromic  point  from  the  Atlantic  Ocean  to  the  Caribbean  Sea 
and  back  to  the  Atlantic,  one  recognizes  a major  tidal  distortion  caused  by  ocean  ridges,  which  has  long  been  dis- 
covered by  practical  tidalists  (see,  e.g.,  Harris.  1904;  Bogdanov,  1961 ; Defant,  1961 ; and  Luther  and  Wunsch,  1974). 
As  the  tide  crosses  the  ridge  between  the  islands,  it  suffers  a distinct  amplitude  jump  and  a significant  phase  shift. 
For  example,  north  of  Puerto  Rico  (PRI)  and  Hispaniola  and  in  the  southeast  around  Barbados  (BB),  the  computed 
and  empirical  Atlantic  tide  data  display  a higher  water  level  and  an  earlier  or,  respectively,  a delayed  cresting  time 
than  the  adjacent  tide  data  on  the  Caribbean  side.  In  particular,  in  full  agreement  with  the  observations,  the  tidal 
retardation  time  can  easily  exceed  30°  (=*  1 hr).  The  distortion  seems  to  depend  on  the  angle  with  which  the  tidal 
crest  spills  over  the  ridge.  Maximum  distortion  appears  to  be  associated  with  a normal  crossing.  It  may  be  pointed 
out  that  the  realistic  resolution  of  tidal  distortions  by  ocean  ridges  (see  below  and  Part  II)  constitutes  probably  the 
most  significant  improvement  of  the  present  model  over  all  earlier  hydrodynamical  models. 

The  Atlantic  portion  of  the  Caribbean-Atlantic  amphidromic  rotation  is  opposed  by  a southward  advancing 
tide  from  about  Newfoundland  in  the  north  and  by  an  eastward  progressing  tide  from  about  Cape  Cod  to  Cape 
Hatteras  in  the  west.  As  a result  of  this  interaction  of  three  opposing  tidal  waves,  the  middle  latitudes  (around 
n = 60°)  of  the  Atlantic  display  very  few  variations  in  tidal  amplitudes  and  phases.  In  the  high-amplitude  sections 
between  Nova  Scotia  and  Cape  Cod  and  between  Cape  Hatteras  and  Florida's  coast,  the  Caribbean-Atlantic  rotation 
system  seems  to  be  less  affected  by  the  opposing  tidal  waves.  The  extreme  tidal  amplitudes  occurring  at  the  former 
shore  section  are  probably  caused  by  the  frontally  advancing  tidal  crest  splashing  against  the  shallow  coastal  corner. 

Although  the  computed  tide  in  the  Gulf  of  St.  Lawrence  displays  the  well-known  amphidromic  point  (Defant, 
1961),  the  grid  system  is  much  too  crude  to  attach  a high  accuracy  to  the  tidal  constants  in  this  border  sea.  For  the 
same  reason,  the  tidal  data  listed  between  Florida,  Cuba,  and  the  Bahamas  are  naturally  less  accurate  than  those  in 
the  open  oceans. 

Tables  6a  and  6b  illustrate  the  smoothness  with  which  the  computed  tide  of  the  northeastern  Pacific  Ocean 
attaches  itself  to  the  empirical  tide  data  along  the  North  American  west  coast.  The  tidal  constants  observed  at  the 
islands  of  Guadalupe  (GI)  and  Farallon  (FI),  at  the  Cobb  Seamount  (CS),  and  at  the  included  near-shore  deep-sea 
stations  fit  realistically  well  into  the  computed  surrounding  tide.  The  amplitudes  and  phases  of  the  excluded  off- 
shore deep-sea  measurements  in  the  Pacific  agree  within  2 cm  and  6°,  respectively,  with  the  computed  data  (Table 
4),  which  is  just  the  same  accuracy  as  in  the  Atlantic. 

Perhaps  the  most  prominent  feature  of  this  area  is  the  amphidromic  point  ®, around  which  the  M2  tide 
rotates  counterclockwise.  This  amphidromic  system  was  predicted  by  Munk  et  al.  (1970)  and  Irish  et  al.  (1971)  in 
almost  identical  geographical  position.  Farlier  hydrodynamical  tide  models  failed  to  resolve  this  system  on  proper 
location,  although  several  models  matched  the  empirical  data  along  the  coast  quite  well.  Since  the  northeastern 
Pacific  falls  short  in  major  bottom  and  coastal  irregularities  when  compared  to  the  northwestern  Atlantic,  the 
indicated  rapid  loss  of  quality  in  westerly  direction  seemed  disappointing.  Yet,  as  will  be  demonstrated  below,  this 
shortcoming  could  have  been  concluded  from  the  obvious  failure  of  those  models  to  reasonably  reproduce  the 
tide  over  most  of  the  north  and  central  Pacific  Ocean. 

As  was  mentioned  before,  the  author's  preliminary  tide  model  (Schwiderski,  1976)  used  a bathymetry  that 
failed  to  represent  the  hydrodynamical  barrier  effects  of  the  Marianas.  Nampo,  Kuril,  Aleutian,  and  Hawaiian 
Ridges,  as  well  as  of  other  seamount  chains.  Consequently,  the  M2  tide  of  almost  the  whole  central,  western,  and 
northern  Pacific  area  was  modeled  as  a single  huge  amphidromic  system,  as  pictured  by  the  similar  maps  of  other 


numerical  tidalists  such  as  Zahel  (1971)  and  Estes  (1975,  1977).  The  clockwise-rotating  Pacific  tide  was  free  to 
sweep  undisturbed  into  the  Philippine,  Okhotsk,  and  Bering  Seas.  By  the  time  the  computed  tidal  crest  reached 
the  Aleutian  Islands,  it  was  just  about  180°  out  of  phase.  When  the  original  bathymetry  was  replaced  by  hydrody- 
namically  defined  depth  data  (Sect.  5.B),  the  entire  Pacific  Ocean  resembled  a whirlpool  after  some  continued 
computations  over  several  quarter  periods.  The  amphidromic  system  weakened  and  its  center  slipped  slowly  south- 
ward, but  drastically  improved  phases  appeared  gradually  along  the  Aleutian  Ridge  confirming  the  anticipated 
effect  of  ocean  ridges. 

The  complete  turnaround  of  the  Pacific  M2  tide  near  the  Aleutian  Islands  was  speeded  up  when  the  empirical 
tidal  constants  were  introduced  into  the  model.  In  fact,  a repeat  of  the  same  computations  settled  the  Pacific  Ocean 
tide  into  its  final  position  in  a rather  dramatic  fashion.  Striking  improvements  were  registered  over  the  whole  Pacific 
and,  of  course,  also  over  the  Atlantic  and  Indian  Oceans. 

As  is  depicted  in  Tables  7a  and  7b  for  the  north-central  Pacific,  the  amphidromic  system  is  replaced  by  a 
low-amplitude  tide.  It  appears  to  be  locked  in  between  the  Aleutian  and  Hawaiian  Ridges  in  the  north  and  south 
and  also  between  the  Emperor  Seamount  chain  in  the  west  and  the  high-amplitude  tide  in  the  east,  which  progresses 
in  a westerly  direction  from  the  west  coast  of  North  America  (Tables  6a  and  6b).  The  amplitude  topography  of 
this  area  resembles  the  low-amplitude  tide  in  the  Caribbean  Sea  (Table  5a).  When  the  westward-advancing  tidal 
wave  enters  the  region  between  the  Aleutian  and  Hawaiian  Ridges,  it  suffers  a remarkable,  almost  symmetric  retarda- 
tion at  both  ridges.  In  fact,  as  the  visible  (0°  = 360°)  cotidal  line  in  Table  7b  reveals,  the  crest  front  of  the  tidal  wave 
assumes  the  shape  of  an  almost  symmetric  wedge.  If  one  traces  the  0°  phase  line  westward  beginning  at  both  ridges, 
one  can  infer  a definite  idea  about  the  realistic  reproduction  of  the  tide  in  this  region.  At  both  ends,  the  0°  phase  is 
in  full  agreement  with  the  empirical  data.  As  the  observed  phases  grow  westward  along  both  ridges,  so  grow  propor- 
tionally the  distances  of  the  0°  phase  line  from  the  ridges. 

The  present  computed  Mj-tide  model  indicates  no  longer  any  symptoms  of  the  original  phase  problems  (see 
above)  at  the  Aleutian  and  Hawaiian  Ridges.  The  computed  amplitudes  and  phases  approach  the  empirical  tidal 
constants  from  both  sides  of  the  ridges  as  smoothly  as  could  be  desired.  As  the  tidal  wave  spills  over  both  ridges 
in  northwestward  or  southwestward  directions,  respectively,  it  suffers  a tidal  distortion  similar  to  that  found  before 
at  the  Caribbean  Ridge.  Amplitude  jumps  and  major  phase  shifts  are  again  in  complete  agreement  with  observations 
(see  the  remarks  of  Luther  and  Wunsch,  1974).  It  is  particularly  gratifying  to  find  the  phase  shift  well  developed 
along  the  whole  length  of  the  Hawaiian  Ridge  from  the  Island  of  Hawaii  to  Midway,  even  though  only  few  empirical 
data  were  used  at  both  ends.  Also,  it  may  be  noticed  that  the  observed  tidal  constants  at  the  distant  and  isolated 
island  stations  of  Pribilof  (PF),  Midway  (MW),  and  Johnston  (JI)  are  all  realistically  well  integrated  by  the  surround- 
ing computed  tide. 

Ironically,  the  old  and  new  Mj-tide  maps  constructed  by  Bogdanov  (1961)  and  Luther  and  Wunsch  (1974) 
by  pure  intuition  and  simple  rules  of  thumb  from  empirical  data  came  closest  to  the  present  charts.  Indeed,  both 
maps  display  no  amphidromic  system  in  the  north-central  Pacific.  As  will  be  verified  in  Part  II,  the  computed  amphi- 
dromic points  between  the  Cook  and  Society  Islands  and  near  the  southern  edge  of  the  Solomon  Islands  are  both  in 
almost  identical  positions  with  those  charted  by  the  same  authors.  Nevertheless,  their  detailed  distribution  of  ampli- 
tudes and  phases  is  still  significantly  different  from  the  present  one. 

Perhaps  the  most  spectacular  display  of  the  high  quality  of  both  the  computed  and  the  observed  tidal  data  is 
brought  out  by  Tables  8a  and  8b  depicting  the  high-amplitude  tide  of  the  central  Pacific.  Indeed,  unlike  any  other 
open  ocean  area,  the  tabulated  region  is  dotted  with  numerous  tide  gauge  stations  at  island  groups  and  at  scattered 
isolated  islands.  In  addition  to  the  fully  listed  island  chains,  there  are  the  isolated  islands:  Johnston  (JI),  Wake  (WI), 
Kusaie  (KI),  Ocean  (Ol),  Funafuti  (FI),  Wallis  (IW),  Niue  (NI),  and  Norfolk  (NF).  The  corresponding  observed 
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Table  6a.  M2  Tidal  Amplitudes,  £ (cm),  of  the  Northeastern  Pacific  Ocean 
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Table  6b.  M2  Tidal  Phases,  6 (°),  of  the  Northeastern  Pacific  Ocean 
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Table  7b.  M2  Tidal  Phases,  6 ( ° ),  of  the  North-Central  Pacific  Ocean 
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Table  8a.  M2  Tidal  Amplitudes,  £ (cm),  of  the  Central  Pacific  Ocean 
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Table  8b.  M2  Tidal  Phases,  6 ( 0 ),  of  the  Central  Pacific  Ocean 
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tidal  constants  listed  in  nongeographical  arrangement  appeared  incoherent  and,  hence,  uncorrelated,  giving  rise 
to  doubt  their  true  value.  Yet,  the  computed  tidal  wave  sweeps  across  the  whole  area  in  a southwesterly  direction 
with  little  variation  of  its  high  amplitude.  As  the  wave  crest  passes  through  the  many  checkpoints,  it  integrates  and 
correlates  without  a single  exception  all  the  empirical  data  into  one  coherent  unity. 


C.  CONCLUSIONS 

The  hydrodynamical  modeling  of  global  ocean  tides,  described  and  tested  above,  shows  that  the  original 
question  posed  by  contemporary  researchers  (Section  l)can  be  successfully  answered.  In  fact,  one  infers  from  the 
evaluation  of  the  constructed  M2  tide  (Sect.  6.B)  that  it  is  now  possible  to  predict  the  M2-tide  elevation  of  the 
ocean  surface  over  the  geoid  anywhere  in  the  open  oceans  with  an  accuracy  of  better  than  5 cm.  This  accuracy 
goes  considerably  beyond  the  originally  desired  error  bounds.  It  leaves  ample  room  for  superposable  errors  due 
to  the  additional  tidal  constituents  listed  in  Table  1,  which  must  be  constructed  with  equivalent  relative  accuracy. 
The  computation  of  the  leading  three  components,  S2,  Kj , and  Op  is  presently  in  progress. 

Naturally,  the  achieved  high  accuracy  of  the  M2  tide  in  the  open  oceans  drops  off  somewhat  near  continental 
or  island  stations  where  empirical  data  are  missing  or  are  less  accurate  themselves  (see  the  introduction  to  the  British 
Admiralty  Tide  Tables,  1977).  Also,  less  accurate  predictions  must  be  anticipated  in  small  border  seas,  bays,  estu- 
aries, and  channels  where  the  1°  by  1°  grid  system  precludes  a sufficient  resolution.  To  improve  the  present  tide 
model  in  those  areas,  significantly  improved  observations  will  be  needed  along  with  a locally  refined  network  and 
corresponding  bathymetric  data. 
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NOAA/Atlantic  Oceanographic  and  Meteorological  Lab. 
Physical  Oceanography  Laboratory 
ATTN:  G.  A.  Maul 
H.  M.  Byrne 


NOAA/Pacific  Marine  Environmental  Lab. 

Seattle,  WA  98105 
ATTN:  Dr.J.  R.  Apel 
H.  0.  Mofjeld 

C.  A.  Pearson 

NOAA/Geophysical  Fluid  Dynamics  Laboratory 
Princeton  University 
Princeton,  NJ  08540 
ATTN:  Dr.  J.  Smagorinsky 
Dr.  K.  Bryan 
Dr.  M.  D.  Cox 

NOAA/National  Center  for  Atmospheric  Research 
Boulder,  CO  80303 
ATTN:  Dr.  W.  R.  Holland 

NASA/Goddard  Space  Flight  Center 
ATTN:  Dr.  J.  W.  Siry 
Dr.  P.  Musen 

D.  E.  Smith 
J.  G.  Marsh 

NASA/Wallops  Station 
Information  Processing  and  Analysis  Branch 
Wallops  Island,  VA  23337 
ATTN:  C.  D.  Leitao 

Director 

U.S.  Army  Ballistic  Research  Laboratory 
Aberdeen  Proving  Ground,  MD  21005 
ATTN:  DRDAR-TSB-S  (STINFO) 

Smithsonian  Astrophysical  Observatory 
60  Garden  St. 

Cambridge,  MA  02138 
ATTN:  Dr.  E.  M.  Gaposchkin 
Dr.  G.  C.  Wiffenbach 
Ms.  Dianne  Hills 

National  Science  Foundation 
1951  Constitution  Ave.,  N.W. 

Washington,  DC  20550 

ATTN:  Mathematical  Sciences  Division 


Scripps  Institution  of  Oceanography 
University  of  California,  San  Diego 
La  Jolla,  CA  92037 
ATTN:  Dr.  W.  H.  Munk 

Dr.  M.  C.  Hendershott 
Prof.  B.  D.  Zetler 
Prof.  S.  M.  Smith 
Prof.  H.  W.  Menard 

MIT/Dept.  Earth  & Planetary  Sciences 
Cambridge  MA  02139 
ATTN:  Dr.  C.  Wunsch 
Dr.  D.  S.  Luther 

Woods  Hole  Oceanographic  Institute 
Woods  Hole,  MA  02543 
ATTN:  Dr.  H.  M.  Stommel 
Dr.  G.  Veronis 
Dr.  N.  P.  Fofonoff 

Battelle  Columbus  Laboratories 
505  King  Ave. 

Columbus,  OH  43201 
ATTN:  A.  G.  Mourad 

Dr.  J.  W.  Chamberlain 

Editor,  Rev.  Geophys.  and  Space  Physics 

Rice  University 

Houston.  TX  77001 

Dr.  R.  H.  Rapp 
Ass.  Editor,  JGR 
Ohio  State  University 
Dept,  of  Geodetic  Science 
1958  Neil  Ave. 

Columbus,  OH  43210 

Dr.  R.  O.  Reid 

Texas  A&M  University 

College  Station,  TX  77843 

Florida  State  University 
Dept,  of  Oceanography 
Tallahassee,  FL  32306 
ATTN:  Dr.  J.  J.  O'Brien 
Dr.  W.  Stur  ,es 


Dr.  R.  H.  Estes 

Business  and  Technological  Systems.  Inc. 
Aerospace  Building.  Suite  605 
10210  Greenbelt  Rd. 

Seabrook.  MD  20801 

T.  V.  Martin 

Sci.  Res.  and  Appl.  Group 

Washington  Analytical  Services  Center.  Inc. 

6801  Kenilworth  Ave. 

Riverdale.  MD  20840 


Dr.  S.  K.  Jordan 

The  Analytic  Sciences  Corporation 
6 Jacob  Way 
Reading.  MA  01867 

The  Rand  Corporation 

Santa  Monica.  CA  90406 

ATTN.  Director.  Climate  Program 

The  Hydrographer  of  the  Navy 
Ministry  of  Defence 
Hydrographic  Department 
Taunton.  Somerset,  England 

Director 

Bureau  Hydrograpltique  International 
7.  Ave.  President  J.  F.  Kennedy 
Monte-Carlo.  Monaco 

G.  C.  Dohler.  P.  Eng. 

Canadian  Hydrographic  Services 
Ocean  and  Aquatic  Affairs 
615  Booth  St. 

Ottava.  Canada  K I A0H3 


Prof.  D.  E.  Cartwright 
Institute  of  Oceanographic  Sciences 
Bidston  Observatory 
Birkenhead.  L437RA,  U.K. 

lnstitut  fur  Meereskunde 
Universitat  Hamburg 
Heimhuder  Str.  71 
2 Hamburg  1 3,  Germany 
ATTN:  Prof.  Dr.  W.  Hansen 
Dr.  W.  Zahel 

Dr.  G.  Seeber 

Tech.  Universitat  Hannover 
Nienburger  Str.  6 
D3  Hannover.  Germany 

Dr.  J.  Schmitt 

Astron.  Geodasie  & Satellitengeodasie 
Techn.  Hochshule  Darmstadt 
Petersenstr.  13 
D-61  Darmstadt.  Germany 

Prof.  P.  Melchior 
Observatoire  Royal  de  Belgique 
Avenue  Circulate  3 
1 180  Brussels.  Belgium 

E.SMETS 

The  Laboratory  of  Hydraulic  Research 

Berchemlei  1 18 

B-2200  Borgerhout.  Belgium 

Dr.  C.  LeProvost 
lnstitut  d.  Mecanique 
BP.  53 

38041  Grenoble.  France 

Meteorological  Res.  Institute 
4-35-8  Koenji-Kita 
Suginami 
Tokyo.  Japan 
ATTN:  T.  Ueno 

M.  Miyazaki 


R.  S.  Chugh.  Director 
Geodetic  Research  Branch 
Dehra  Dun,  India 


Prof.  G.  Obenson 
University  of  Lagos 
Lagos,  Nigeria 


Local:  D 


K,  KOI , K02,  K04,  K05H,  K05J 

K10,  K10LK11.K12,  K12H,  K12S,  K12W 

K20,  K30,  K55,  K60,  K70,  K73,  K73B,  K73S 

R04 

E41 

X210  (2) 

X210  (G1DEP)  (2) 

KOSS  (100) 


1 


[ 
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INFORMATION 


LIST  OF  CORRECTIONS  TO  NSWC/DL  TR-3866 


Page  3: 

Page  22: 
Page  27: 

Page  29: 
Page  31: 
Page  32: 


In  Table  1,  change  anplitude  K of  S2  to:  K(S2>  = 0.112841 
and  D to  : D = d + 365 (y  - 1975)  + Int  [(y-1973)/4] 

In  Eq.  14  change  nQ  to:  nQ  = 1sK(3cos20  - 2)  oos  (at  + }C 
Top  line  change  6 to:  B = 0.90 

In  first  equation  of  Eq.  32  add:  Hx=  Hx/H,  Hxx  = Hxx/H. 

2 

Change  Eq.  37c  to  : *sK(3  sin  3-2). 

In  V-equation  of  Eq.  (43)  change  sin  9 to  sin  26. 

Change  Eq.  45  to: 


n = *sK(3sin  6 - 2)  e 
= K. (3sin29  - 2)  e 
U = bK2  sin  6 e 10t 

V ="*^!^aR  sin  26e'''ot 


iat 

iot 


In  Equation  45a  change  k2  and  to: 

k2  = -^/^3 

K,  = -3aGH  K Jj6BGH  - a2R2  + fiRK,)  + io(6A„  + B R2/H  )]. 

1 O O O O O OJ 

Note:  A higher  order  approximation  of  improved  accuracy  is  in  preparation. 


Page  33: 
Page  47: 

Page  52: 
Page  56: 


Change  to:  = -3K. 

In  B change  for  v=  0: 

m,n  ^ 

-3  sin  4nA6  to  +6sin4nA6. 

In  Eg.  84  change  the  left-side  C to  c. 

In  Eq.  94  change  k.  and  k_  to:  k.  = .06,  k2  = 


.03. 


